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Abstract 


The  vibration  suppression  effectiveness  of  a  flexible  in-line  marine  cable  vibration 
absorber  is  studied.  The  transfer  matrix  method  is  used  to  build  various  numerical  models 
of  vibration  absorbers  in  marine  cable  systems.  The  models  determine  cable  system 
natural  frequencies,  mode  shapes  and  modal  damping  ratios. 

The  introduction  of  absorber  damping  is  shown  to  result  in  complex  roots  to  the 
modal  characteristic  equations.  A  computer  complex  root  solver  is  used  to  solve  for  the 
complex  roots  of  the  characteristic  equations,  resulting  in  complex  system  natural 
frequencies.  The  significance  of  complex  natural  frequencies  is  explained.   Complex 
natural  frequencies  are  used  to  calculate  modal  damping  ratios. 

The  models  demonstrate  that  absorber  effectiveness  is  heavily  dependent  on 
absorber  location,  absorber  mass  and  absorber  length.  Parametric  variation  is  used  to 
achieve  maximum  effectiveness  of  the  flexible  in-line  absorber.  Even  under  optimum 
conditions,  it  is  shown  that  the  absorber  provides  insufficient  damping  to  reduce  vortex- 
induced  vibrations  in  water. 

The  same  transfer  matrix  method  is  used  to  evaluate  the  effectiveness  of  a  mass- 
spring-dashpot  type  absorber  in  a  marine  cable  system.  This  type  of  absorber  is  shown  to 
produce  adequate  damping  to  reduce  vortex-induced  vibrations  in  water.   The  transfer 
matrix  method  used  in  this  thesis  is  validated  by  analyzing  the  same  system  using  an 
approach  by  Den  Hartog  [1].  The  transfer  matrix  approach  combined  with  complex  root 
solving  capability  is  shown  to  provide  an  effective  analysis  method  for  marine  cable 
systems. 
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Chapter  1 
Introduction 

Long  slender  marine  structures  such  as  cables  and  pipelines  are  subject  to  a 
condition  called  vortex  induced  vibration.  This  occurs  when  the  cylindrical  structure  is 
subjected  to  a  cross  flow  current  as  shown  in  figure  1-1: 

Flow  speed,  U  (m/s) 


Vortices  in  wake 

Figure  1-1:  Vortex  Shedding 

At  a  critical  value  of  flow  speed,  U,  a  periodic  shedding  of  vortices  occurs  on  the  upper 
and  lower  surfaces  of  the  cylinder.  This  behavior  is  known  as  a  Von  Karman  street.  As 
each  vortex  is  shed  into  the  wake,  a  localized  pressure  drop  occurs.  The  result  is  a 
harmonically  varying  vertical  (lift)  force  on  the  cylinder.  The  frequency  at  which  vortex 
shedding  occurs  is  related  to  cylinder  diameter  and  flow  velocity  by  a  dimensionless 
parameter  called  the  Strouhal  number: 

f  D 

St=       TJ  (i-i) 

where:  f  =  vortex  shedding  frequency  (Hz) 

D  =  cylinder  diameter  (m) 
U  =  flow  speed  (m/sec) 


The  Strouhal  number  is  approximately  0.2  for  a  wide  range  of  Reynolds  numbers. 
Vortex-induced  vibrations  are  especially  severe  when  the  vortex  shedding  frequency  is 
near  a  resonance  frequency  of  the  structure.  Cross-flow  displacement  of  the  cylinder  at 
this  "lock-in"  frequency  may  be  on  the  order  of  one  cylinder  diameter.   Such  vibrations 
can  fatigue  marine  structures  and  cause  eventual  failure. 

This  problem  motivates  the  investigation  into  a  dynamic  absorber  for  long 
cylindrical  marine  structures.  Before  analyzing  the  effect  of  various  damping  devices,  it  is 
worth  while  to  determine  the  magnitude  of  damping  which  must  be  introduced  into  a 
marine  structure  to  cause  significant  reduction  in  vortex  induced  vibrations.   In  figure  1-2 
below,  Griffin  [13]  shows  how  cross  flow  displacement  varies  with  a  dimensionless 
parameter  called  the  reduced  damping,  SG.  Reduced  damping  is  the  ratio  of  damping 
force  to  fluid  exciting  force. 
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Figure  1-2:  Vibration  Displacement  vs.  Reduced  Damping 


where:  Sq  =  8  71"  St  ^  (  ms  /  Pf  D")  =  reduced  damping  (1-2) 

St  =  Strouhal  number 

ms  =  linear  density  of  cable  (kg/m) 

Pf=  fluid  density  (kg/m") 

D  =  cable  diameter  (m) 
The  reason  displacement  amplitude  approaches  a  plateau  value  at  low  reduced  damping 
(values  below  Sg  about  0.2)  is  due  to  the  self-limiting  nature  of  vortex  induced  vibrations. 
At  values  of  displacement  amplitude  greater  than  approximately  one  cylinder  diameter, 
the  extreme  motion  of  the  cylinder  begins  to  disturb  the  formation  of  the  Von  Karman 
street.   Stuctural  damping  is  not  the  limiting  factor  at  low  values  of  reduced  damping.   The 
reduced  damping  is  valuable  in  determining  the  amount  of  damping  which  must  be  added 
to  the  structure  to  have  an  effect  on  displacement  amplitude.   It  would  not  be  worthwhile 
to  add  damping  such  that  reduced  damping  would  not  exceed  the  "knee  in  the  curve" 
value  of  0.2  illustrated  in  figure  1-2.  This  is  also  a  very  convenient  approach  since  the 
damping  ratio  used  to  calculate  reduced  damping  is  the  damping  which  one  would 
measure  in  an  in  vacuo  transient  decay  test.  Thus,  various  absorber  configurations  can  be 
evaluated  without  the  need  to  consider  hydrodynamic  effects. 

This  thesis  will  study  various  vibration  absorber  configurations  for  marine  cables. 
Most  of  the  analysis  will  focus  on  the  type  of  absorber  shown  in  figure  1-3  on  the 
following  page: 


cable 


Rigid  masses 

Figure  1-3:  Cable  Vibration  Absorber 
This  absorber  was  first  proposed  by  Vandiver  and  Li  [10]  .  Energy  dissipation  occurs  by 

way  of  relative  angular  motion  at  the  pivot  between  the  two  rigid  masses.  The  damping 

device  is  integral  to  the  pivot  but  it  is  shown  schematically  in  figure  1-3.  Butler  [7] 

showed  this  type  of  damping  device  to  be  possible  for  marine  cable  applications.   It  was 

seen  as  a  very  versatile  configuration  for  marine  cables  since  such  an  absorber  could  be 

clamped  on  at  a  desired  location  while  deploying  a  cable.  Or,  if  the  absorber  were 

permanently  installed,  its  size  and  orientation  might  permit  the  absorber  to  be  reeled  onto 

a  cable  drum  when  the  cable  is  not  deployed. 

Chapter  2  describes  the  transfer  matrix  approach  for  analyzing  multi-component 
mechanical  systems  with  finite  boundaries.  It  is  a  very  convenient  method  for  studying 
the  dynamic  behavior  of  marine  cable  systems  with  attached  damping  devices.  The 
analysis  utilizes  the  absorber  transfer  matrix  derived  by  Li  [5]. 

In  chapter  3,  the  transfer  matrix  method  is  used  to  determine  natural  frequencies, 
mode  shapes  and  damping  ratios  of  the  cable  /  absorber  system.  A  method  is  formulated 
for  determining  the  optimum  value  of  damping  constant. 

In  chapter  4  ,  system  parameters  such  as  absorber  length  and  mass  are  varied  to 
study  their  effect  on  damping  performance.   Also,  use  of  multiple  absorbers  and  the  effect 
of  absorber  location  relative  to  vibration  nodes  and  antinodes  was  studied.   Use  of  a  more 


conventional  absorber  type  is  studied  by  merging  the  transfer  matrix  approach  with  a  more 
classical  approach  (Den  Hartog  [1])  with  some  interesting  results. 

Chapter  5  concludes  the  study  and  suggests  areas  for  further  research. 
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Chapter  2 

The  Transfer  Matrix  Method 

The  transfer  matrix  method  allows  for  a  convenient  matrix  representation  of  the 
marine  cable  system.  In  1963,  Pestel  and  Leckie  [2]  recognized  that  the  introduction  of 
numerical  computing  made  matrix  methods  a  powerful  tool  for  studying  mechanical 
system  dynamics.  In  marine  cable  systems,  features  such  as  rigid  bodies  and  discrete 
damping  devices  are  well  suited  to  matrix  modeling  techniques.  Today's  computing  speed 
allows  complicated  marine  cable  systems  to  be  modeled  with  relative  ease. 

2.1  The  Marine  Cable  System 

2.1.1  Assumptions 

It  is  important  here  to  state  the  assumptions  employed  in  representing  the  cable 
system  in  this  way.  These  assumptions  are  those  of  the  classic  linear  string  wave  equation. 
They  were  used  by  Li  [5]  in  developing  the  transfer  matrices  used  in  this  analysis: 

1)  the  cable  systems  are  taut  and  undergo  only  small  deflections 

2)  cable  bending  stiffness  is  neglected 

3)  the  masses  in  the  cable  system  are  rigid 

4)  the  weight  of  these  masses  is  small  compared  to  cable  tension 

5)  the  dynamic  variation  in  tension  in  the  cable  is  small  and  neglected 

6)  only  planar  motion  is  considered  (2  dimensional) 
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2.1.2  Transfer  Matrix  Representation  of  a  Marine  Cable  System 


Transfer  matrices  are  used  to  relate  the  physical  state  of  one  location  in  a 
mechanical  system  to  that  at  another  location  in  the  system.  The  state  at  any  location  in 
the  system  is  represented  by  the  state  vector.  In  a  cable  system,  the  state  at  any  location 
can  be  represented  by  the  state  vector,  z: 


[z] 


'y(x) 

TT(x) 


(2-1) 


where  y  is  the  Fourier  transform  of  cable  transverse  displacement  and  y  is  the  Fourier 
transform  of  cable  slope.  The  product  of  cable  tension  T  and  y  is  a  value  equal  to  the 
transverse  force  caused  by  tension  at  the  location  of  the  state  vector.  Thus,  the  cable 
system  state  vector  contains  a  displacement  value  and  a  force  value.   If  assumption  2  in 
section  2.1.1  were  not  valid  and  cable  bending  stiffness  could  not  be  neglected,  the  state 
vector  would  require  four  elements:  displacement,  slope,  shear  force  and  bending  moment. 
The  use  of  Fourier  transforms  is  useful  since  the  types  of  motion  studied  here  will  limited 
to  free  and  forced  harmonic  vibrations. 

Within  a  mechanical  system,  the  state  vector  at  any  location  can  be  obtained  by 
multiplying  the  state  vector  at  some  other  location  in  the  system  by  the  transfer  matrix  of 
the  component(s)  between  the  two  locations.  In  a  marine  cable  system,  the  transfer 
matrices  are  2  by  2  matrices  represented  as  follows: 


[T]  = 


til    tl2 

t:i    t:2 


(2-2) 
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The  elements  of  the  transfer  matrices  are  derived  using  the  equations  of  motion.   Since 
they  relate  the  Fourier  transforms  of  displacement  and  force,  the  transfer  matrix  elements 
are  functions  of  circular  frequency,  co.  The  transfer  matrix  of  a  component  is  general  in 
that  it  does  not  change  when  the  component  is  installed  in  other  systems.  Cable  system 
components  modeled  here  using  transfer  matrices  include  cable  segments,  lumped  and 
distributed  masses  and  vibration  absorbers. 

Figure  2-1  illustrates  how  transfer  matrices  are  used  to  relate  the  state  vectors  at 
individual  locations  in  the  cable  system: 


/ 

[Zl] 


[Ti] 
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absorber 
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[2j  =  [Tl]  [Zl] 

[z3]  =  [Ti]  [T2]  [zj 

[z?]  =  [T2]  [z2] 

[z(x)]  =  [T1(x)][z,];x<x: 


(2-3) 
(2-4) 
(2-5) 
(2-6) 


Figure  2-1:  Transfer  Matrix  Representation  of  Cable  /  Absorber  System 

The  transfer  matrix  for  the  uniform  cable  section,  [TJ,  comes  directly  from  the 
solution  to  the  linear  string  wave  equation  and  is  given  by: 


[Ti]  = 


cos(k/)  l/kTsin(k/) 

-kTsin(k/)    cos(k/) 


(2-7) 


where:  k  =  wave  number  =  co  /  c  (rrf1) 

c  =  wave  propagation  speed  =  (T  /  p)  "    (m/s) 

T  =  cable  tension  (N) 

p  =  cable  linear  mass  density  (kg/m) 

/  =  cable  length  from  xi  to  x2  (m) 
For  locations  on  the  cable  to  the  left  of  x2,  substitution  of  x  for  /  in  equation  2-7  will  yield 
the  transfer  matrix  [Ti(x)]  as  in  equation  2-6. 

The  transfer  matrix  for  the  absorber,  [T2],  is  given  in  Appendix  B  and  in  the  thesis 
by  Li  [5].  Appendix  B  also  contains  the  transfer  matrices  for  various  vibration  absorbers 
and  discrete  masses  which  may  be  located  within  a  cable  system. 
2.1.3  The  Overall  Transfer  Matrix 

The  product  of  the  three  transfer  matrices,  [Ti]  [T2]  [T3],  equals  the  overall 
transfer  matrix  of  the  system.  For  the  system  shown  in  figure  2-1,  the  boundaries  are 
walls.  For  a  wall,  cable  displacement  would  be  zero  and  cable  force  would  be  unknown. 
Thus,  the  state  vectors  at  the  boundaries  can  be  represented  by: 


W  = 


0 
Ty'(x) 


(2-8) 


The  overall  system  transfer  matrix  can  be  used  to  relate  the  two  boundary  state  vectors: 
[Toverall]  =  [TJ  [T2]  [T3]  (2-9) 

[Z4]  =  [Toven.1,]  [Zi]  (2-10) 
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0 

Ty'M 


tn    ti: 
t:i   t:: 


0 

TT(x.) 


(2-11) 


The  zero  displacement  at  the  boundaries  results  in: 
(tn  xO)+(ti:xTy(xi))  =  0 


(2-12) 


^>ti:=0  (2-13) 

This  is  the  characteristic  equation  for  the  system.   Since  transfer  matrix  elements  are 
functions  of  frequency,  the  roots  of  ti2  of  [T0veraii]  are  the  natural  frequencies  of  the 
system. 
The  mode  shapes  of  the  system  can  then  be  plotted  by  setting  a  unity  force  at  the  wall: 


H  = 


(2-14) 


It  follows  that: 


[z(x)]  =  [Ti(x)]  [zj  ;x<x2 

=  [Ti][T2][T3(x)][Zl];x>X3 


(2-15) 
(2-16) 


Note  that  [T3(x)]  is  expressed  by  substituting  the  quantity  (x  -  x? )  for  /  into  equation  2-7. 
This  coordinate  transformation  is  required  because  transfer  matrix  expressions  normally 
reference  a  coordinate  system  which  is  located  on  the  component  that  the  transfer  matrix 
represents. 
In  general, 


[z(x)]  =  [T(x)]  [Zl]  - 


tn(x)  ti:(x) 
t:.(x)  t::(x) 


(2-17) 


where  [T(x)]  is  the  transfer  matrix  product  of  the  system  components  between  the  left  wall 
and  x. 
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The  first  element  of  [z(x)j  is  therefore: 

y(x)  =  t.:(x)  (2-18) 

Equation  2-18  specifies  the  displacement  mode  shape  of  the  system.  Note  that  ti2(x)  in 
equation  2-18  is  an  element  of  the  transfer  matrix  between  the  left  wall  and  an  arbitrary 
point,  x,  and  is  not  the  same  as  ti:  of  T0veraii  in  equation  2-13. 

In  chapter  3,  this  approach  will  be  applied  to  a  specific  marine  cable  /  absorber 
system. 
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Chapter  3 

Vibration  Absorber  Analysis 

As  mentioned  in  the  introduction,  the  vibration  absorber  which  will  be  the  primary 
focus  of  this  analysis  is  that  shown  below  in  figure  3-1 : 


<*— >< c ^ *- 


h  =  c  +  d 


l  =  e  +  f 


Figure  3-1:  In  Line  Vibration  Absorber 

The  transfer  matrix  expression  for  this  absorber  was  derived  by  Li  [5]  and  is  quite  lengthy, 

requiring  a  page  and  a  half  to  state  (see  Appendix  B).  This  is  because  the  transfer  matrix 
must  account  for  both  the  translational  and  rotational  motion  of  the  masses,  damping  and 
stiffness.  MATLAB  was  used  to  compute  the  values  of  this  transfer  matrix  and  was  used 
for  all  other  modeling  and  calculation  in  this  thesis.  Appendix  A  contains  the  MATLAB 
programs  STRTNG1  and  STRING2  which  calculate  the  transfer  matrices  of  the  absorber 
and  cable  segments,  respectively. 

3.1  Modeling  the  Vibration  Absorber  in  a  Cable  Svstem 

Since  the  nature  of  this  analysis  was  numerical,  a  baseline  absorber  system  of  the 
type  shown  in  figure  3-1  was  modeled.  System  parameters  such  as  cable  length,  cable 
diameter  and  absorber  mass  were  set  to  values  which  could  realistically  be  expected  in  a 
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marine  system.  Optimum  absorber  performance  can  be  obtained  by  systematically  varying 
these  parameters.  Baseline  system  parameters  are  as  follows: 

Total  system  length:  50.8  m 

Cable  tension,  T:  5000  N 

Cable  density,  p:  1  kg/m 

Absorber  link  mass,  M:  9  kg 

Absorber  link  moment  of  inertia,  I:      0. 12  kg  m2 

a  and  b  dimensions:  0.03  m 

c,  d,  e,  f  dimensions:  0.2  m 

Total  absorber  length,  Lat,s:  0.8  m 

Stiffness,  K:  lON/m 

Damping,  R:  varies  N/(m/s) 

Note  that  the  stiffness,  K,  is  relatively  insignificant.  Cable  tension  is  the  dominant  factor  in 
providing  a  restoring  force  for  the  absorber. 

Absorber  analysis  consisted  of  the  following  steps: 

1)  Obtain  natural  frequencies 

2)  Determine  damping  ratio 

3)  Determine  displacement  mode  shapes 

4)  Vary  system  parameters  to  optimize  damping  ratio 
This  process  is  detailed  in  the  following  sections. 


18 


3.2  The  Cable  /  Absorber  System  With  no  Damping 

It  is  valuable  to  first  analyze  the  system  and  demonstrate  the  analysis  process  with 
the  damping  value  of  the  absorber  set  to  zero.  This  keeps  all  calculations  in  the  real 
domain  and  provides  for  a  better  understanding  of  the  introduction  of  damping. 

3.2.1  System  Natural  Frequencies 

Chapter  2  detailed  the  process  for  determining  the  characteristic  equation  and 
natural  frequencies  using  transfer  matrices.  The  overall  transfer  matrix  of  the  system, 
[Toveraii],  is  computed  as  follows: 

[Toverall]  =  [Ti]    [T2]    [T?]  (3-1) 

where: 

[Ti]  and  [T:,]  are  the  transfer  matrices  of  the  left  and  right  cable  segments 
as  in  equation  2-7 

[T2]  is  the  transfer  matrix  of  the  absorber  as  in  Appendix  B 
Due  to  the  fixed  wall  boundary  conditions,  the  system  characteristic  equation  is 

ti2=0  (3-2) 

where  ti2  is  the  first  row,  second  column  element  of  T0Veraii 

A  MATLAB  script  file,   STRTNG4,  was  used  to  generate  figure  3-2  on  the 
following  page,  a  plot  of  ti2  versus  frequency,  co.  The  zeroes  oft,:  are  the  natural 
frequencies  of  the  system. 
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T(1.2)  CHARACTERISTIC  FOR  SYSTEM  WITH  ZERO  DAMPING 


Figure  3-2:  ti2  Characteristic  of  Undamped  System 

3.2.2  System  Mode  Shapes 

Once  the  system  natural  frequencies  are  known,  transfer  matrix  multiplication  can 
be  used  to  determine  the  displacement  and  force  mode  shapes.  This  is  accomplished  by  a 
MATLAB  program  called  STRTNG6  which  takes  the  system  frequency  as  an  input.  The 
state  vector  at  the  left  wall  is  set  as  follows 


Zicft  — 


(3-3) 


where  the  first  element  indicates  the  zero  displacement  at  the  wall  and  the  second  element 
is  set  at  one.  This  results  in  a  mode  shape  which  is  normalized  to  a  unity  transverse  force 
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at  the  system  boundaries.  The  displacement  mode  shape  is  calculated  along  the  length  of 
the  system  as  described  in  Chapter  2: 
from  left  wall  to  absorber: 


[z(x)]  =  [T,(x)][z,ett]  = 


y(x)  =  tp.(x), 


tn    ti: 
t:i    t:: 


0 


y(x) 

Ty'(x) 


where  ti2(x)  is  an  element  of  [Ti(x)] 
from  absorber  to  right  wall: 
y(x)  =  ti:(x) , 


(3-4) 


(3-5) 


(3-6) 


where  t^(x)  is  an  element  of  the  matrix  product  [Ti]  [T2]  [T?(x)]. 
Similarly,  the  transverse  force  mode  shape  would  equal  the  t22  element  of  the  transfer 
matrix  as  a  function  of  distance  from  the  left  wall.  The  resulting  displacement  mode 
shapes  for  the  first  three  vibration  modes  are  shown  in  figure  3-3  on  the  following  page. 
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Figure  3-3:  Undamped  System  Displacement  Mode  Shapes 
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3.3  Introduction  of  Absorber  Damping 

The  introduction  of  absorber  damping  takes  the  transfer  matrix  calculations  into 
the  complex  domain.  In  addition  to  natural  frequency  and  mode  shape,  damping  ratio  can 
now  be  determined.  Damping  ratio  is  the  parameter  used  in  determining  absorber 
effectiveness.  Damping  ratio  is  used  to  compute  the  reduced  damping  in  figure  1-2  and  to 
determine  if  sufficient  damping  exists  to  reduce  vortex-induced  vibrations. 
3.3.1  System  Natural  Frequencies 

The  same  approach  which  was  used  to  compute  undamped  natural  frequencies  in 
section  3.2. 1  is  used  here  to  compute  damped  natural  frequencies.   As  shown  in  section 
3.2. 1,  the  characteristic  equation  of  the  system  is: 

ti2=0  (3-7) 

where  tn  is  an  element  of  the  overall  system  transfer  matrix.  Figure  3-4  on  the  following 
page  shows  the  magnitude  of  ti2  as  a  function  of  frequency,  co.  Note  that  when  damping 
is  present,  Xn  is  a  complex  quantity.  The  figure  also  shows  ti2  for  the  undamped  system 
in  the  same  frequency  range 

The  damped  curve  is  quite  different  in  appearance  from  the  undamped  case.  For 
alternating  modes,  the  damped  value  of  ti2  fails  to  reach  a  zero  for  real  values  of  natural 
frequencies.  This  behavior  is  significant  and  will  later  be  shown  to  have  a  very  physical 
explanation. 
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Figure  3-4:   Magnitude(t12)  Versus  Frequency  for  Damped  vs.  Undamped  System 
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In  order  to  obtain  all  the  zeroes  of  ti2,  complex  values  of  frequency  must  be  considered. 
This  requires  the  use  of  a  computer  complex  root  solver,  MULLER  (Appendix  A), 
developed  by  Conte  and  de  Boor  [14]  and  converted  into  MATLAB  format  by  Wilson  and 
Roberts.  Before  proceeding  further  with  the  results  of  the  complex  root  search,  the 
meaning  of  a  complex  natural  frequency  must  first  be  discussed. 

The  harmonic  motion  at  a  point  in  the  cable  system  can  be  described  by: 
y(x)  =  Ae-  (3.8) 

where  A  is  displacement  amplitude. 
For  imaginary  natural  frequencies: 

co  =  coreai  +  ico,mag  (3-9) 

Substituting: 

y(x)  =  AeH",real  +  K!,imag)t  (3-10) 

—  ^  gi(o>real)t  e-(wimag)t  (3-11) 

=  A[cos(oW)  +  isin(corealt)]  e<miwa**       (3-12) 

The  real  part  of  which  is: 

=  A  cos(cDreait)  e*°imas)t      (3-13) 
So,  we  see  that  for  free  vibrations  the  damped  system  vibrates  at  coreai  with  a  decay 
envelope  e^fflUMg)t.  The  ratio  of  the  real  and  imaginary  parts  of  the  natural  frequency  is  the 
modal  damping  ratio: 

C  =  cOimag/coreai  (3-14) 

Now  that  the  significance  of  a  complex  natural  frequency  is  understood,  we  can 
return  to  the  analysis  of  the  cable  /  absorber  system.  The  complex  root  solver,  MULLER, 
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is  iterative  and  requires  a  starting  guess  for  the  complex  frequency  root.  The  undamped 
natural  frequencies  were  used  as  initial  guesses  when  running  MULLER.  With  the 
assistance  of  Dr.  Rama  Rao,  MIT,  a  routine  called  MULRUN  (Appendix  A)  was 
developed  which  successively  runs  MULLER  for  many  values  of  damping  constant  to 
obtain  an  output  which  demonstrates  the  optimum  value  of  damping.  Figure  3-5  below 
shows  the  output  of  MULRUN  for  the  first  vibration  mode  of  the  baseline  system. 
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Figure  3-5:  Demonstration  of  Optimum  Damping  Value 
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The  MULRUN  output  shows  that  as  damping  constant,  R,  is  varied,  there  is  an 

intermediate  damping  value  which  optimizes  (maximizes)  damping  ratio.  It  was  theorized 

that  the  reason  damping  ratio  decreases  at  high  values  of  damping  constant  is  that  relative 

motion  between  the  absorber  masses  decreases  above  the  optimum  value.  This  was 

confirmed  by  creating  a  transfer  matrix  model  of  a  uniform  bar  (see  Appendix  B)  in  place 

of  the  vibration  absorber  in  the  baseline  system  as  shown  in  figure  3-6: 

damper 


Rigid  masses 


I   I 


Uniform  bar 
cable 


Figure  3-6:  Representing  High  Stiffness  Absorber  as  Uniform  Bar 

The  uniform  bar  had  the  same  dimensions  and  mass  as  those  of  the  absorber.   It  was 
shown  that  the  cable  system  with  the  uniform  bar  at  its  center  had  the  same  natural 
frequency  and  mode  shape  as  the  absorber  with  very  high  damping  constant.   So  at  very 
high  damping,  relative  rotational  motion  between  the  absorber  links  is  restricted  and  the 
absorber  is  unable  to  dissipate  energy.  The  absorber  acts  as  a  uniform  bar. 
3.3.2  System  Mode  Shapes  With  Damping 

Once  the  system  natural  frequencies  are  known,  transfer  matrix  multiplication  can 
be  used  to  determine  the  displacement  and  force  mode  shapes  using  the  same  approach  as 
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for  the  undamped  case.  This  is  accomplished  by  a  MATLAB  program  called  STRJNG6 
which  takes  system  frequency  as  an  input.  For  the  damped  case,  complex  natural 
frequencies  are  input  into  STRJNG6.  The  state  vector  at  the  left  wall  is  again: 


[zieft] 


(3-15) 


The  displacement  mode  shape  is  calculated  along  the  length  of  the  system  as  in  section 
3.2.2.    The  resulting  mode  shapes  for  the  first  three  vibration  modes  are  shown  in  figure 
3-7. 
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Figure  3-7:  Damped  System  Displacement  Mode  Shapes 
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Inspection  of  these  mode  shapes  reveals  that  for  the  even  mode,  the  absorber  is  located  at 
a  modal  antinode.  There  is  no  relative  motion  between  the  absorber  links.  The  absorber 
simply  rotates  about  its  center  as  if  it  were  a  uniform  bar.  This  is  confirmed  further  by 
using  MULRUN  for  the  same  modes  to  determine  optimum  damping  ratio.  In  figure  3-8 
below,  it  is  seen  that  for  the  even  vibration  mode,  damping  ratio  is  insignificant  even 
though  damping  is  present.  Relative  motion  between  the  absorber  links  is  required  if  the 
absorber  is  to  dissipate  energy. 
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Fi«ure  3-8:  Damping  Effect  in  First  Three  Vibration  Modes 
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So,  from  initial  analysis  the  following  observations  are  made: 

1)  An  optimum  value  of  damping  constant  is  one  which  adds  sufficient 
damping  while  not  overly  restricting  motion  between  the  two  absorber 
links. 

2)  Absorber  effectiveness  is  very  dependent  on  absorber  location. 
Absorbers  located  at  modal  antinodes  provide  the  most  energy  dissipation. 

These  observations  are  consistent  with  the  behavior  of  more  conventional  types  of 
vibration  absorbers,  such  as  those  installed  in  discrete  (vice  continuous)  mechanical 
systems. 

In  chapter  4,  parameters  such  as  absorber  mass,  absorber  length  and  number  of 
absorbers  will  be  varied  to  study  their  effect  on  absorber  performance. 
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Chapter  4 

Parameters  Affecting  Absorber  Performance 

In  previous  chapters,  it  has  been  demonstrated  how  natural  frequencies,  mode 
shapes  and  damping  ratios  can  be  obtained  for  a  vibration  absorber  in  a  marine  cable 
system.  It  was  shown  how  the  damping  constant,  R,  can  be  varied  to  arrive  at  an 
optimum  damping  ratio.  The  tools  are  now  in  place  to  study  the  effects  of  other 
parametric  variations  on  absorber  performance.  Improvements  in  absorber  performance 
will  be  sought  by  varying  the  following: 

1)  Absorber  length 

2)  Absorber  mass 

3)  Number  of  absorbers  installed  in  a  single  cable  system 

4)  Absorber  location 

4.1   Required  Damping 

Before  proceeding  with  further  attempts  to  optimize  damping,  it  would  be  helpful 
to  know  what  value  of  damping  is  required  in  the  system  to  reduce  vortex  induced 
vibrations.  Recall  from  figure  1-2  that  there  is  a  value  of  reduced  damping  below  which 
structural  damping  provides  no  assistance  in  reducing  vortex  induced  vibrations. 
Structures  with  reduced  damping  in  this  plateau  range  will  have  displacement  amplitudes 
which  are  restricted  by  limit-cycle  reductions  to  the  lift  forces.   Thus,  there  is  no  gain  in 
adding  small  amounts  of  damping  which  will  not  raise  reduced  damping  above  that  value 
required  to  reduce  the  response  to  less  than  that  resulting  from  limit  cycle  behavior.  From 
figure  1-2,  this  threshold  value  of  reduced  damping  is  as  follows: 
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So  >  0.2 
where: 

Sg  =  8tT  St2C(ms/pfD2)  (4-1) 

For  the  baseline  system  in  water  (see  section  3.1),  this  results  in  a  threshold  (required) 
value  of  damping  ratio: 

C  >  04  (4-2) 

The  highest  value  of  damping  ratio  achieved  thus  far  (see  figure  3-8)  is  approximately 
Q  =  .0012.   So,  we  see  that  if  the  absorber  being  studied  is  to  be  of  use  in  reducing  vortex 
induced  vibrations  in  water,  further  optimization  is  required. 

4.2  Variation  of  System  Parameters 

4.2.1  Absorber  Length 

In  figure  3-7,  mode  shapes  for  the  first  three  modes  were  shown.  These  low 
modes  have  very  long  wavelengths,  resulting  in  the  situation  shown  in  figure  4- 1(a)  on  the 
following  page.  When  absorber  length  is  small  compared  to  wavelength,  relative  angular 
motion  between  the  absorber  links  is  very  small,  even  in  odd  numbered  vibration  modes. 
Since  energy  dissipation  relies  on  this  relative  motion,  short  absorber  links  result  in  poor 
absorber  performance.   As  absorber  length  is  increased,  relative  motion  between  the 
absorber  links  is  increased,  and  energy  dissipation  increases  as  shown  in  figure  4- 1(b).   A 
key  dimensionless  parameter  of  use  here  would  be  the  ratio  of  absorber  length  to  the  half 
wave  length  of  a  given  mode.  This  parameter  is  designated  as  the  length  ratio: 
LR  =  2Labsorber/A.  (4-3) 
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Figure  4-1:  Effect  of  Absorber  Length  on  Damper  Motion 

Variation  of  absorber  length  was  accomplished  in  an  indirect  way.  Instead  of  varying  the 

absorber  link  length,  the  effect  of  absorber  length  was  studied  by  varying  X,  or  simply  by 
comparing  the  damping  ratios  of  different  modes  of  the  baseline  system.   Starting  with 
mode  1,  optimum  damping  values  were  determined  for  each  mode.   As  discussed  in 
section  3.3.2,  even  numbered  vibration  modes  exhibited  no  damping  effects.  The  best 
absorber  performance  (as  determined  by  maximum  damping  ratio)  was  achieved  in  the 
29th  mode. 
In  this  mode: 

?i  =  3.50m 

Labsorbei  =-8  m 

LR=457 

In  modes  beyond  mode  29,  higher  length  ratio  resulted  in  more  restricted  absorber  motion 
and  reduced  damping  ratios. 


j  j 


Figure  4-2  shows  the  damping  value  optimization  and  mode  shape  for  mode  29. 
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Figure  4-2:   Damping  Optimization  and  Mode  Shape  for  Mode  29 

It  can  be  seen  that  with  an  optimum  value  of  damping  constant,  damping  ratio  is 

approximately  1.7%.  This  falls  well  short  of  the  value  of  4%  required  to  reduce  vortex 
induced  vibration  amplitude  in  water. 
4.2.2  Absorber  Mass 

In  further  efforts  to  increase  the  amount  of  damping  which  the  absorber  adds  to 
the  system,  the  mass  of  the  absorber  links  was  varied.  A  useful  dimensionless  parameter 
for  this  analysis  was  the  mass  ratio: 

MR  =  mabsorber  /  (Labsorber  pcable  )  (4-4) 

This  is  the  ratio  of  the  absorber  mass  to  the  mass  of  the  cable  segment  displaced  by  the 
absorber. 
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The  variation  of  mass  ratio  was  performed  for  the  29th  mode  of  vibration  of  the 
cable  system.  This  mode  was  selected  because  it  this  was  the  mode  for  which  optimum 
absorber  length  was  shown  .  Thus  the  result  of  this  variation  represents  an  attempted 
optimization  of  both  absorber  mass  and  length.  There  is  no  guarantee  that  optimum  length 
ratio  remains  at  0.457  when  mass  ratio  is  varied.  Figure  4-3  shows  how  optimum 
damping  ratio  varies  with  absorber  mass: 
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Figure  4-3:  Affect  of  Mass  Ratio  on  Mode  29  Damping 


Figure  4-3  shows  that  increasing  absorber  mass  ratio  tends  to  increase  optimum 
damping  ratio.  Also,  lower  values  of  damping  constant  are  required  to  achieve  the 
optimum  damping  value  when  a  higher  mass  ratio  is  used.  This  figure  demonstrates  that 
the  baseline  system  (MR=22.5)  has  the  optimum  mass  ratio  because  it  results  in  the 
highest  system  damping  ratio.  When  mass  ratio  was  increased  beyond  22.5,  the  system 
shifted  to  an  even  numbered  vibration  mode  where  damping  is  insignificant. 

If  this  type  of  analysis  were  to  be  performed  on  a  real  marine  system,  one  would 
have  to  exercise  care  in  changing  absorber  mass.  Changing  absorber  mass  changes  the 
mass  of  the  vibrating  system  and  thus  changes  system  natural  frequencies.  The  marine 
system  would  be  designed  with  some  range  of  current  velocity  and  corresponding  vortex 
shedding  frequencies  in  mind.  If  one  were  not  careful,  one  could  change  system  mass  such 
that  the  absorber  was  optimized  at  some  frequency  outside  the  range  of  concern. 
4.2.3  Behavior  at  Realistic  Values  of  Absorber  Mass 

Admittedly,  a  mass  ratio  of  MR=22.5  corresponds  to  a  very  large  absorber  mass, 
equaling  36%  of  the  total  cable  mass.  Absorbers  of  this  extreme  size  were  considered 
only  in  an  attempt  to  achieve  the  required  damping  level  of  section  4. 1.  A  more  typical 
absorber  mass  might  be  no  more  than  10%  of  the  cable  mass.  This  would  correspond  to  a 
mass  ratio  of  MR=6.25.  In  this  section,  an  absorber  of  this  size  was  analyzed  over  many 
vibration  modes  and  corresponding  length  ratios.  The  results  are  shown  in  figure  4-4  and 
4-5  on  the  following  page. 
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OPTIMIZED  DAMPING  RATIO  VS.  MODE 
NUMBER  AND  LENGTH  RATIO 
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Figure  4-4:  Optimized  Damping  Ratio  vs.  Mode  Number  and  Length  Ratio  for  Mass  Ratio  =  6.25 
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Figure  4-5:   Optimized  Damping  Constant  vs.  Mode  Number  and 
Length  Ratio  for  Mass  Ratio  =  6.25 

Figures  4-4  shows  that  an  optimum  damping  ratio  is  achieved  in  mode  53  with  a  length 
ratio  of  LR=0.83.  Figure  4-5  shows  how  required  damping  decreases  dramatically  with 


increasing  length  ratio.  This  behavior  was  discussed  in  section  4.2.1 


4.3  Number  of  Absorbers 

Up  to  this  point,  variation  of  system  parameters  has  not  resulted  in  achieving  the 
required  damping  ratio  of  Section  4.1.  It  would  be  very  convenient  if  one  could 
multiplicatively  increase  damping  ratio  by  simply  increasing  the  number  of  absorbers 
installed. 

This  addition  of  absorbers  was  performed  for  mode  3.  A  transfer  matrix  model 
was  constructed  for  a  cable  system  with  vibration  absorbers  at  its  mode  3  antinodes.   An 
absorber  mass  ratio  of  MR=6.25  was  used.   All  other  system  parameters  were  those  of  the 
section  3.1  baseline  system.  The  MATLAB  programs  which  model  this  system  are  located 
in  directory  MULTABS  in  Appendix  A    These  programs  are  functionally  identical  to 
those  used  to  evaluate  the  system  with  one  absorber. 

Figure  4-4  on  the  following  page  shows  the  mode  3  mode  shapes  for  the  system 
with  one  absorber  and  the  system  with  three  absorbers.  For  both  cases,  a  mass  ratio  of 
MR=6.25  and  a  length  ratio  of  LR=.047  was  used.  The  maximum  damping  value 
achieved  in  each  case  is  indicated. 
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Figure  4-6:  Comparison  of  Single  Absorber  to  Multiple  Absorbers 

Inspection  of  figure  4-6  reveals  that  optimized  damping  ratio  roughly  doubled 
when  the  number  of  absorbers  was  tripled.   Thus,  damping  ratio  is  not  multiplicative  with 
the  number  of  absorbers  installed  in  the  system.   One  should  also  observe  from  fisure  4-4 
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that  mode  3  natural  frequency  decreases  considerably  with  the  additional  mass  loading  of 
the  extra  absorbers. 

4.4  Absorber  Location  for  Operation  in  Multiple  Modes 

In  previous  sections,  absorber  performance  has  been  looked  at  one  mode  at  a  time. 
Parameters  such  as  damping  constant  and  absorber  mass  have  been  varied  and  the  effect 
on  modal  damping  ratio  has  been  observed.  This  has  improved  understanding  of  the 
system  but  it  is  not  the  way  in  which  a  marine  vibration  absorber  system  would  be 
designed.   In  a  system  subject  to  vortex  induced  vibrations,  some  range  of  potential 
currents  would  exist.  This  range  of  currents  would  be  directly  correlated  to  an  excitation 
bandwidth  by  the  Strouhal  relationship  of  equation  1-1.   So,  an  absorber  design  problem 
would  be  concerned  with  multiple  modes.  This  section  will  illustrate  how  the  proper 
location  of  a  single  absorber  can  optimize  performance  in  multiple  vibration  modes. 

In  section  3.3.2,  it  was  observed  that  the  most  effective  location  for  an  absorber 
was  at  a  modal  antinode.  For  mode  3,  these  locations  would  be  L/6,  L/2,  and  5L/6.  Now 
suppose  we  wanted  to  design  an  absorber  system  which  would  be  effective  for  modes  3 
and  4.  The  mode  4  antinodes  are  at  L/8,  3L/8,  5L/8.   The  two  modes  share  no  common 
antinodes.  In  fact,  the  L/2  antinode  for  mode  3  is  a  mode  4  node.   So,  an  absorber  located 
at  L/4  offers  no  absorber  effect  for  mode  4. 

A  compromise  was  arrived  at  by  constructing  a  model  of  a  system  with  an 
absorber  at  7L/16.  This  location  is  halfway  between  the  3L/8  mode  4  antinode  and  the 
L/2  mode  3  antinode.  The  absorber  used  had  a  mass  ratio  of  MR=6.25.  The  MATLAB 
programs  for  modeling  and  analyzing  this  system  are  located  in  directory  OFFCENTER  in 
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Appendix  A.  Table  4-1  below  compares  absorber  effectiveness  for  the  compromise 
absorber  location  versus  absorber  locations  optimized  for  a  single  mode: 


Absorber  Located  for  Mode  3 
Effectiveness 

Absorber  Located  for  Dual 
Mode  Effectiveness 

Mode  3  Damping 
Ratio 

.00028 

.00023 

Mode  4  Damping 
Ratio 

0 

.00023 

Table  4-1:  Dual  Mode  Absorber  Optimization 

The  results  indicate  that  relocation  of  the  absorber  causes  some  loss  of  mode  3 
damping  ratio  while  raising  mode  4  damping  ratio  from  zero  to  a  value  equal  to  that  of 
mode  3.  This  is  the  type  of  tradeoff  would  have  to  be  performed  in  the  design  of  a  real 
marine  cable  absorber  system. 

Figure  4-7  on  the  following  page  shows  the  mode  3  and  4  mode  shapes  with  the 
absorber  located  at  7L/16.  This  figure  demonstrates  the  utility  of  the  STRTNG6  mode 
plotter  in  visualizing  unorthodox  mode  shapes. 

So,  absorbers  can  be  located  for  effectiveness  in  both  even  and  odd  numbered 
modes.  In  choosing  absorber  location,  the  designer  of  a  marine  cable  system  would  have 
to  consider  the  excitation  bandwidth  and  vibration  modes  within  that  bandwidth.  If  a 
single  absorber  location  could  not  provide  energy  dissipation  in  all  possible  excitation 
modes,  then  multiple  absorbers  would  need  to  be  installed. 
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Figure  4-7:  Mode  3  and  4  Shapes  with  Absorber  at  7L/16 
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4.5  Consideration  of  a  Different  Absorber  Type 

Up  to  this  point  in  the  analysis,  considerable  insight  has  been  gained  in  the 
behavior  of  the  hinged,  in-line,  vibration  absorber  in  the  marine  cable  system.  However, 
this  type  of  absorber  has  not  been  able  to  produce  the  necessary  damping  ratio  required  to 
reduce  vortex  induced  vibrations  in  a  cable  system  in  water.  (See  section  4. 1 ). 

In  this  section,  a  simpler  and  more  conventional  absorber  will  be  analyzed.  This 
absorber  is  shown  in  figure  4-8  below: 


cable 


S] 


>rin«,  k 


n 


Damper,  R 


Mass,  m 


Fij^jre  4-8:  Han«»in«-Mass-Sprin«  Dashpot  Absorber 

The  absorber  is  of  the  conventional  mass-spring-dashpot  type.  It  is  usually  associated 
with  discrete  mechanical  systems  but  could  also  be  used  to  dissipate  energy  in  a  cable 
system.  Li  and  Vandiver  [10]  derived  the  following  transfer  matrix  for  this  absorber: 


[T]  = 


1 

arm  (jft>R  +  k) 
k  -  co'm  +  jcoR 


(4-5) 


A  numerical  model  of  a  cable  system  using  this  absorber  was  constructed  in 
MATLAB.  The  cable  itself  is  identical  to  that  of  the  baseline  system  of  section  3.1. 
However,  the  absorber  of  figure  4-8  is  located  at  the  center  of  the  cable  in  place  of  the  in- 
line bending  type  absorber.  The  programs  comprising  this  model  are  located  in  directory 
HANG  in  Appendix  A. 

Den  Hartog  [1]  formulated  a  method  for  choosing  optimum  absorber  mass  and 
stiffness  for  this  type  of  absorber.  His  formulation  was  for  discrete  systems  but  modal 
analysis  of  the  continuous  cable  system  allows  Den  Hartog's  method  to  be  applied  in  an 
approximate  way.  The  approach  uses  the  following  parameters: 

u.  =  m  /  M  =  mass  ratio  =  absorber  mass/  main  mass 


CO 


=  V k  /  m  =  natural  frequency  of  absorber 


Q  =  VK  /  M  =  natural  frequency  of  main  system  (cable  without  absorber) 

f  =  co  /  Q  =  frequency  ratio 

Den  Hartog  determined  that  the  correct  tuning  of  the  mass  spring  dashpot 
absorber  results  from  the  frequency  ratio  being  set  as  follows: 

f=l/(l+u)  (4-6) 

This  equation  was  used  to  set  absorber  stiffness    To  do  this,  an  example  mass  ratio  of 
u.  =  0.2  was  selected.  The  main  mass  was  set  at  the  modal  mass  of  the  cable  which  equals 
one  half  of  total  cable  mass.  Q.  was  set  at  the  cable  first  mode  natural  frequency. 
Equation  4-6  then  allowed  computation  of  the  absorber  stiffness,  k 
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The  following  values  were  used  in  the  analysis  of  the  vibration  absorber  for  the  first 
vibration  mode  of  the  cable  system: 

m  =  5  kg  =  absorber  mass 

M  =  25  kg  =  cable  modal  mass  =  Vi  of  total  cable  mass 

u  =  m/M  =  0.2 

Q  =  4.44  rad/sec  =  cable  first  mode  natural  frequency  without  absorber 
Using  equation  4-6: 

f=  0.833  =co/Q 

co  =  3.70  rad/sec  = 


k  =  68.5  N/m,  the  correct  stiffness  for  optimum  tuning 

Once  stiffness  and  mass  were  set,  MULRUN  was  used  to  determine  the  optimum 
damping  constant  for  the  absorber  designed  in  the  procedure  above.  Figure 
4-9  below  shows  the  resulting  damping  ratio  performance  and  first  mode  shape: 


MODAL  DAMPING  RATIO  VS.  DAMPING  CONSTANT 


20  25  30  35 

x  location  (m) 


Figure  4-9:   Dampei  Optimization  for  First  Mode  with  Hanging  Mass  Spring  Dashpot  Absorber 
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It  can  be  seen  that  optimum  damping  ratio  is  approximately  0.22,  far  higher  than 
any  damping  ratio  achieved  with  the  in  line  absorber.  This  damping  ratio  is  also  above  the 
critical  value  of  C,  -  .04  necessary  to  reduce  vortex  induced  vibration  amplitude  in  water. 
This  result  also  provides  a  reasonable  reconciliation  of  results  between  the  Den  Hartog 
method  and  the  transfer  matrix  method  used  in  this  thesis.  Den  Hartog  states  that  the 
optimum  damping  constant  for  the  mass  spring  dashpot  absorber  can  be  determined  by  the 
following  formula: 

— —  =  [Su/CSO+u)3)]1-  (4-7) 

2  Ll  m 

For  a  mass  ratio  of  u.  =  0.2,  this  results  in  an  RoptimumOf  9.25  N/(m/s).  This  corresponds 

roughly  with  the  optimum  damping  constant  of  15  N/(m/s)  shown  in  figure  4-7. 

The  hanging  mass-spring-dashpot  absorber  represents  a  large  improvement  in 

absorber  effectiveness  over  the  bending  type  absorber  discussed  elsewhere  in  this  thesis. 

This  absorber  has  some  possibility  for  use  in  marine  cable  systems  and  further  analysis 

should  be  conducted. 
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Chapter  5 
Conclusion 

In  this  thesis,  the  transfer  matrix  method  was  used  to  construct  various  numerical 
models  of  marine  cable  systems.  The  performance  of  various  vibration  absorber 
configurations  was  studied.  The  following  are  the  main  findings  of  this  research: 

1)  Introduction  of  damping  in  the  transfer  matrix  model  causes  the  system 
characteristic  equation  to  have  complex  roots.  The  use  of  a  computer  complex 
root  solver  was  introduced  and  was  shown  to  be  useful  in  calculating  system 
natural  frequencies  and  damping  ratios. 

2)  Upon  understanding  the  significance  of  complex  natural  frequencies,  computer 
programs  using  transfer  matrices  were  created  to  plot  system  mode  shapes. 
The  ability  to  visualize  mode  shapes  was  very  important  in  understanding 
system  behavior. 

3)  The  use  of  an  in-line  bending  type  vibration  absorber  was  studied  in  some 
depth.  A  method  was  devised  for  tuning  such  an  absorber  to  achieve 
maximum  damping  within  a  marine  cable  system.  Variation  of  parameters  such 
as  absorber  mass,  length  and  number  of  absorbers  was  performed  to  study  the 
effect  of  these  variations  on  absorber  performance.  Even  after  optimization, 
the  in-line  absorber  was  unable  to  demonstrate  sufficient  damping  to  reduce 
vortex-induced  vibrations  in  water. 
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4)   The  use  of  a  hanging  mass-spring-dashpot  absorber  was  also  studied  for  use  in 
a  marine  cable  system.  This  absorber  was  shown  to  have  the  ability  to  create 
sufficient  damping  to  reduce  vortex  induced  vibrations  in  water. 

The  main  focus  of  this  thesis  was  to  study  the  performance  of  the  in-line  bending 
type  absorber  in  a  finite  length  marine  cable  system.   Although  it  was  shown  that  this 
absorber  provides  insignificant  damping  to  reduce  vortex-induced  vibrations  in  water, 
considerable  insight  was  gained  in  this  thesis.  A  working  model  was  created  which  has  the 
ability  to  calculate  system  natural  frequencies,  mode  shapes  and  damping  ratios. 
Additionally,  it  should  be  noted  that  the  in-line  absorber  could  be  useful  in  reducing  vortex 
induced  vibrations  in  air.  The  studies  in  this  thesis  would  be  useful  in  optimizing  the 
absorber  for  this  application. 

Areas  for  potential  research  related  to  this  thesis  include: 

1)  Use  of  the  transfer  matrix  models  in  this  thesis  to  study  other  marine  cable 
system  characteristics,  such  as  the  effects  of  distributed  masses  within  the 
system.  Li  [5]  performed  considerable  work  in  this  area  but  the  ability  to 
visualize  mode  shapes  would  provide  a  good  supplement  his  work. 

2)  Extension  of  the  analysis  procedure  in  this  thesis  to  other  systems  which  can  be 
represented  by  transfer  matrices. 

3)  Further  study  of  the  mass-spring-daspot  absorber  for  use  in  marine  cable 
systems. 
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Appendix  A:  MATLAB  Files 


This  appendix  contains  all  MATLAB  files  used  in  this  thesis    The  files  are  divided 
into  four  directories,  each  corresponding  to  a  different  cable  system  model.    The 
following  is  a  list  of  these  directories  and  a  description  of  the  system  they  model: 

CENTER:  A  cable  system  50  meters  long  with  an  in  line  bending  type  absorber  at 
its  center. 

OFFCENTER:  A  cable  system  50  meters  long  with  an  in  line  bending  type 
absorber  located  at  7L/16. 

MULTABS:  A  cable  system  50  meters  long  with  three  in  line  bending  type 
absorbers,  each  located  at  a  mode  3  antinode. 

HANG:  A  cable  system  50  meters  long  with  a  hanging  mass  spring  dashpot 
absorber  located  at  its  center. 
Variables  and  Dimensions 

The  following  variables  are  used  throughout  the  MATLAB  programs  in  this 
Appendix: 


Variable 

Corresponding  Parameter 

Units 

a,  b,  c,  d,  e,  f 

cable  dimensions 

m 

M1,M2 

absorber  link  masses 

kg 

T,  H 

cable  tension 

N 

K 

stiffness 

N/m 

R 

damping  constant 

N/(m/s) 

11,12 

absorber  link  moments  of  inertia 

kgm2 

P 

cable  mass  density 

kg/m 

c 

cable  wave  propagation 

speed 

m/s 

k 

cable  wave  number 

m"1 

w 

circular  frequency 

rad/s 

50 


CENTER 

The  following  is  a  typical  sequence  for  using  the  model  CENTER  to  determine 
system  natural  frequencies  and  mode  shapes  and  optimize  system  damping: 

1)  In  MATLAB,  change  directory  to  CENTER. 

2)  Edit  STRING  1  for  desired  absorber  properties  and  STRTNG2  for  desired  cable 
properties. 

3)  In  STRTNG4  and  STRING1,  ensure  damping  constant  is  set  to  zero.  Run  STRTNG4 
to  plot  ti2  of  [Toveran]  versus  frequency.  The  zeroes  of  this  plot  are  the  undamped  system 
natural  frequencies. 

4)  If  desired,  plot  undamped  system  mode  shapes  as  follows: 

a)  For  mode  of  interest,  enter  natural  frequency  in  MATLAB  window. 

b)  Run  STRTNG6  to  plot  mode  shape. 

c)  Repeat  for  other  modes  of  interest. 

5)  Perform  absorber  optimization  for  individual  modes  as  follows: 

a)  Edit  MULRUN  to  iterate  over  desired  range  of  damping  constants. 

b)  Also  in  MULRUN,  set  zrs  to  the  undamped  natural  frequency  for  the  mode  of 
interest.  This  serves  as  a  starting  guess  for  the  complex  root  solver. 

6)  Set  damping  constant  to  global  in  STRING1.  Run  MULRUN  to  plot  real  and 
imaginary  natural  frequency  and  damping  ratio  versus  damping  constant.  The  optimum 
damping  constant  is  that  which  maximizes  damping  ratio. 

7)  Plot  damped  system  mode  shape  as  follows: 

a)  Set  damping  constant  to  the  optimum  value  in  STRING  1. 

b)  Enter  the  complex  natural  frequency  in  the  MATLAB  window. 

c)  Run  STRTNG6  to  plot  mode  shape. 

8)  Repeat  steps  5  through  7  for  other  modes  of  interest. 
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%  STRING  1 

%  This  file  contains  the  transfer  matrix  for  a  vibration 
%  absorber. 

%  Initial  data 

T=5000; 

a=.03; 

b=.03; 

c=.2; 

d=2; 

e=2; 

f=.2; 

Ml=9.0; 

M2=M1; 

H=T; 

Il=.12; 

12=11; 

K=10; 

R=4000; 

h=c+d; 

l=e+f; 

al=(T*d)+(H*c)+(K*aA2)-(W2*Il)+(i*w*R*bA2); 

a2=(T*f)+(H*e)+(K*a"2)-(w^2*I2)+(i*w*R*b^2); 
a3=(K*a~2)+(i*w*R*b~2); 

gl=(-(Ml*l+M2*f)*al- 
(Ml*c*a3)+(w'2*Ml*M2*f*c"2))/((M2*e*al)+(Ml*d+M2*h)*a3+(w"2*Ml*M2*c*e*d)); 

g2=(-(M2*f*h'2)- 
(Ml*l*d^2)+(al*l+a3*h)/w^2)/((M2*e*al)+(Ml*d+M2*h):f:a3+(w'2*Ml*M2*c*e*d)); 

((Ml*l*d*c)+(al*l+a3*h)/w~2)/((M2*e*al)+(Ml*d+M2*h)*a3+(w*2*Ml*M2*c*e*d)); 

g4=-((M  1  *1+M2*f)*a3+(M  1  *c*a2)+(w~2*M  1  *M2*e*c*f))/((a2*h+a3*l)/wA2- 
(Ml*d*r2+M2*r2*h)); 

g5=((a2*h+a3*l)/w^2+(M2*e*f*h))/((a2;l:h+a3:,:l)/w'2-(Ml*d*r2+M2:i:r2*h)); 

g6=(-a2:,:(Ml*d+M2*h)-(M2:l:e*a3)+(w^2*Ml*M2*e'2*d))/((a2*h+a3*l)/w^2- 
(Ml*d*r2+M2*r2*h)); 


52 


t(l,l)=gl+g3*((g4+gl*g6)/(l-g3*g6)); 
t(l,2)=g2+g3*((g5+g2*g6)/(l-g3*g6)); 
t(2,l)=(g4+gl*g6)/(l-g3*g6); 
t(2,2)=(g5+g2*g6)/(l-g3*g6); 


D3 


%  STRING2 


%  This  file  calculates  the  transfer  matrix  for  a  string  of  length  L. 


%Initial  data 

T=5000; 

L=25; 

p=1.0; 

c=sqrt(T/p); 

k=w/c; 

t(l,l)=cos(k*L); 

t(l,2)=sin(k*L)/(k*T); 

t(2,l)=-k*T*sin(k*L); 

t(2,2)=cos(k*L); 
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%  STRING4 

%  This  file  plots  the  overall  T(l,2)  for  a  system  with  3  components. 

W=0; 

Ttotl2=0; 

R=0; 

for  n=  1:250; 

w=.l*n; 

W(n)=w; 

string2 

Tl=t; 

T3=t; 

string  1 

T2=t; 

Ttot=Tl*T2*T3; 

Ttotl2(n)=Ttot(l,2); 

end 

plot(W,abs(Ttotl2)) 

grid 

xlabel('w(rad/sec)') 

ylabel('T(l,2)') 


%  STRING5 

%  This  file  contains  the  transfer  matrix  for  a  uniform  bar. 
%    Initial  data: 

%  bar  length 
lb=8; 

M=18.0; 

T=5000; 
Ib=M*lb*w~2/T; 

k=l/(2*(6+Ib)); 

t(l,l)=k*4*(3-Ib); 

t(l,2)=k*12*lb/T; 

t(2,l)=-k*M*wA2*(12-Ib); 

t(2,2)=k*4*(3-Ib); 
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%  STRING6 

%         THIS  FILE  PLOTS  THE  MODE  SHAPES  OF  A  STRING 

%         WITH  A  DAMPER  AT  ITS  CENTER.  BEFORE  RUNNING, 

%         ENTER  FREQUENCY  IN  THE  MATLAB  WINDOW  AND 

%         ENTER  DAMPING  IN  STRING  1 . 

y=0; 

f=0; 

T=5000; 

p=1.0; 

c=sqrt(T/p); 

k=w/c; 

x=0:.  1:50.8; 

for  n=  1:251 

tl(l,l)=cos(k*x(n)); 

tl(l,2)=sin(k*x(n))/(k*T); 

tl(2,l)=-k*T*sin(k*x(n)); 

tl(2,2)=cos(k*x(n)); 

y(n)=tl(l,2); 

f(n)=tl(2,2); 

end 


string  1 

ts(l,l)=cos(25*k); 

ts(l,2)=sin(25*k)/(k*T); 

ts(2,2)=cos(25*k); 

ts(2,l)=-k*T*sin(25*k); 

forn=258:509 

tr(l,l)=cos(k*(x(n)-25.8)); 

tr(  1 ,2)=sin(k*(x(n)-25.8))/(k*T); 

tr(2,l)=-k*T*sin(k*(x(n)-25.8)); 

tr(2,2)=cos(k*(x(n)-25.8)); 

tf=ts*t*tr; 

y(n)=tf(l,2); 

f(n)=tf(2,2); 

end 

plot(x(  1 :25 1  ),real(y(  1 :25 1  )),x(258:509),real(y(258:509))) 

grid 

xlabel('x  location') 
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ylabel('real(displacement  mode  shape)') 
title('3rd  MODE  MODE  SHAPE') 


58 


%  STRING  14 

%  This  file  computes  T(l,2)  for  a  system  with  3  components. 

function  T12=stringl4(w) 

string2 

Tl=t; 

T3=t; 

string  1 

T2=t; 

Ttot=Tl*T2*T3; 

T12=Ttot(l,2); 
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%  MULLER 

function  [zros]=muller(fn,nprev,np,maxit,epl,ep2,fnreal,zros); 
%  zeros=muller(fn,nprev,maxit,ep  1  ,ep2,fnreal,zeros) 
%  This  function  is  called  by  the  driver  program  mulrun.m 
%  Determines  up  to  np  zeros  of  the  function  specified  by  fn, 
%  using  quadratic  interpolation,  i.e.,  Muller's  Method. 
%  This  function  is  a  translation  of  the  Fortran  routine  given 
%  by  S.  D.  Conte  and  Carl  de  Boor,  "Elementary  Numerical 
%  Analysis",  McGraw  Hill  Book  Company,  1980.  The  conversion 
%  was  made  by  Howard  Wilson  and  Chris  Roberts,  Engineering 
%  Mechanics  Department,  University  of  Alabama. 
% 

%  To  find  more  than  one  zero,  function  muller  uses  a  procedure 
%  known  as  deflation,  which  is  implemented  in  function  dflate. 
%  fn:  fuction  that  we  want  to  find  the  roots. 

%  zros:  an  array  of  initial  estimates  of  all  desired  roots,  set  to  zero  if 
%        no  better  estimates  are  available. 
%  epl:  relative  error  tolerance  on  the  root. 
%  ep2:  error  tolerance  on  the  function. 
%  maxit:  maximum  number  of  iterations  per  root  allowed. 
%  np:  number  of  roots  desired. 

%  nprev:  number  of  roots  previously  computed,  normally  zero. 
%  fnreal:  a  logical  number;  if  TRUE,  the  program  forces  all  approximations  to 
%  all  the  roots  to  be  real.  This  makes  it  possible  to  use  this  routine 

%  even  if  f(x)  is  defined  only  for  real  x. 

%  Here  we  can  take  fnreal=l,  if  we  only  want  to  fine  the  real  roots;  we 

%  can  take  fnreal=0  if  we  want  to  find  both  the  real  and  complex  roots. 

fn='stringl4'; 

eps  1  =max(ep  1 , 1  .e- 12);  eps2=max(ep2, 1  .e-20); 
for  i=nprev+l:np;  kount=0; 
while  1 

%    Compute  the  first  three  estimates  for  zero  as 
%    zero,  zero+0.5,  zero-0.5 

zero=zros(i);  h=0.00001+j*0.00001;  zp=zero+h;  zm=zero-h; 

[fzr,dvdflp,kount,zrosl=dflate(fn,zp,i,kount,zros); 

[fzr,fzrprv,kount,zros]=dflate(fn,zm,i,kount,zros); 

hprev=-2.0*h;  dvdflp=(fzrprv-dvdflp)/hprev; 

[fzr,fzrdfl,kount,zros]=dflate(fn,zero,i,kount,zros); 

while  2 

%      Do  while  kount<maxit  or  h  is  relatively  big  or  fzr=f(zero) 
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%      is  not  small  or  fzrdfl=fdeflated(zero)  is  not  small  or  not 
%      much  bigger  than  its  previous  value  fzrprv. 

di vdf  1 =(fzrdfl-fzrprv)/h;  di vdf2=(divdf  1  -dvdf  1  p)/(h+hprev) ; 
hprev=h;  dvdflp=divdfl;  c=divdfl+h*divdf2; 
sqr=c*c-4.0*fzrdfi*divdf2; 

if  fnreal*real(sqr)<0.0,  sqr=0.0;  else,  sqr=sqrt(sqr);  end; 
if  real(c)*real(sqr)+imag(c)*imag(sqr)<0,  den=c-sqr; 
else,  den=c+sqr;  end 
if  abs(den)<0.0,  den=1.0;  end 
h=-2.0*fzrdfl/den;  fzrprv=fzrdfl;  zero=zero+h; 
if  kount>maxit,  break,  else 
while  3 
[fzr,fzrdfl,kount,zros]=dflate(fn,zero,i,kount,zros); 
if  zros(i)==zero*  1.001,  break;  end 
%  Check  for  convergence 

if  abs(h)<epsl*abs(zero),  break;  end 
if  max(abs(fzr),abs(fzrdfl))<eps2,  break;  end 
%  Check  for  divergence 

if  abs(fzrdfi)>10.0*abs(fzrprv),  h=h/2;  zero=zero-h; 
else,  break;  end 
end;  %  while  3 
end;  %  end  kount  <=  maxit 
if  zros(i)==zero*  1.001,  break;  end 
if  abs(h)<epsl*abs(zero),  break;  end 
if  max(abs(fzr),abs(fzrdfl))<eps2,  break;  end 
end;  %  while  2 
if  kount>maxit,  break;  end 
if  abs(h)<epsl*abs(zero),  break;  end 
if  max(abs(fzr),abs(fzrdfl))<eps2,  break;  end 
end;  %  while  1 
zros(i)=zero; 
end;  %  for  loop 
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%  MULRUN 


C7     ************************************** 

x=0; 

y=0; 

global  R 

%  input  specs  for  Muller 

for  n=  1:51 

R=100*(n-1); 

x(n)=R; 

%    tolerance,  name  of  function  that  computes  T12  etc: 

zrs  =121.;  %  Starting  guess  for  the  root,  can  be  real  or  complex 

fn  =  'string  14';  %  name  of  program  that  computes  T12 
maxit  =  500;  %  max  #  of  iterations 

epl       =       le-5;  %tolerance       on       the       root       (start       with       lower 

%  values,  le-2or  so.) 

tolfcn  =  le-5;  %  relative  tolerance  on  T12 

fnreal  =0;  %  =  0  complex  search,  =  1,  real  search 

nprev  =  0;  %  leave  these  as  0 

np  =  1 ;  %  and  1 


ep2  =  tolfcn*stringl4(zrs);  %  This  sets  the  tolerance  on  Tl 2  to  be 

%       le-5      of      T12      estimated      at      your      starting 
%  guess 

%  Call  Muller 

root  =  muller(fn,nprev,np,maxit,epl,ep2,fnreal,zrs); 
y(n)=root; 
end 
zeta=imag(y)./real(y); 

subplot(311) 

plot(x,imag(y)) 

grid 

title('4th  mode') 

xlabel('R') 

ylabel('imag(nat  freq)') 

subplot(312) 

plot(x,real(y)) 
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grid 
xlabel('R') 

ylabel('real(nat  freq)') 

subplot(313) 

plot(x,zeta) 

grid 

xlabel('R  (N/(m/s))') 

ylabelC modal  damping  ratio') 


%  DFLATE 

%  This  is  the  subroutine  needed  by  muller. 

function  [fzero,fzrdfl,kount,zros]=dflate(F,zero,i,kount,zros); 

%  (*fzero,fzrdfl,kount,zeros*)=dflate(F,zero,i,kount,zeros) 

%  This  function  which  is  called  by  function  muller  is  translated 

%  from  the  Fortran  routine  given  by  S.  D.  Conte  and  Carl  de  Boor, 

%  "Elementary  Numerical  Analysis",  McGraw  Hill  Book  Company,  1980. 

%  The  translation  was  made  by  Howard  Wilson  and  Chris  Roberts, 

%  Engineering  Mechanics  Department,  University  of  Alabama. 

kount=kount+l;  fzero=feval(F,zero);  fzrdfl=fzero; 

if  i<2,  return;  end 

forj=2:i 

den=zero-zros(j-l); 

if  abs(den)==0.0 
zros(i)=zero*  1.001;  return 

else 
fzrdfl=fzrdfl/den; 

end 
end 


64 


OFFCENTER 

The  following  is  a  typical  sequence  for  using  the  model  OFFCENTER  to 
determine  system  natural  frequencies  and  mode  shapes  and  optimize  system  damping: 

1)  In  MATLAB,  change  directory  to  OFFCENTER. 

2)  Edit  STRING  1 A  for  desired  absorber  properties  and  STRJNG2AL  and  STRING2AR 
for  desired  cable  properties. 

3)  In  STRING4A  and  STRING  1  A,  ensure  damping  constant  is  set  to  zero.  Run 
STRING4A  to  plot  ti2  of  [Toveraii]  versus  frequency.  The  zeroes  of  this  plot  are  the 
undamped  system  natural  frequencies. 

4)  If  desired,  plot  undamped  system  mode  shapes  as  follows: 

a)  For  mode  of  interest,  enter  natural  frequency  in  MATLAB  window. 

b)  Run  STRTNG6A  to  plot  mode  shape. 

c)  Repeat  for  other  modes  of  interest. 

5)  Perform  absorber  optimization  for  individual  modes  as  follows: 

a)  Edit  MULRUN  to  iterate  over  desired  range  of  damping  constants. 

b)  Also  in  MULRUN,  set  zrs  to  the  undamped  natural  frequency  for  the  mode  of 
interest.  This  serves  as  a  starting  guess  for  the  complex  root  solver. 

6)  Set  damping  constant  to  global  in  STRING  1  A.  Run  MULRUN  to  plot  real  and 
imaginary  natural  frequency  and  damping  ratio  versus  damping  constant.  The  optimum 
damping  constant  is  that  which  maximizes  damping  ratio. 

7)  Plot  damped  system  mode  shape  as  follows: 

a)  Set  damping  constant  to  the  optimum  value  in  STRING  1  A. 

b)  Enter  the  complex  natural  frequency  in  the  MATLAB  window. 

c)  Run  STRING6A  to  plot  mode  shape. 

8)  Repeat  steps  5  through  7  for  other  modes  of  interest. 
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%  STRING  1 A 

%  This  file  contains  the  transfer  matrix  for  a  vibration 
%  absorber. 

%  Initial  data 

T=5000; 

a=03; 

b=.03; 

c=2; 

d=2; 

e=2; 

f=.2; 

Ml=9; 

M2=M1; 

H=5000; 

Il=.12; 

12=11; 

K=10; 

R=  100000; 

h=c+d; 

l=e+f; 

al=(T*d)+(H*c)+(K*a"2)-(wA2*Il)+(i*w*R*b'2); 

a2=(T*f)+(H*e)+(K*a^2)-(w^2*I2)+(i*w*R*b^2); 
a3=(K*a~2)+(i*w*R*b~2); 

gl=(-(Ml*l+M2*f)*al- 
(Ml*c*a3)+(w'2*Ml*M2*f*c^2))/((M2*e*al)+(Ml*d+M2*h)*a3+(w^2*Ml*M2*c*e*d)); 

g2=(-(M2*f*h^2)- 
(Ml*l*dA2)+(al*l+a3*h)/wA2)/((M2*e*al)+(Ml*d+M2*h)*a3+(wA2*Ml*M2*c*e*d)); 

g3=- 

((M 1  *l*d*c)+(al  *l+a3*h)/w~2)/((M2*e*al)+(Ml  *d+M2*h)*a3+(w~2*M  1  *M2*c*e*d)); 

g4=-((Ml*l+M2*f)*a3+(Ml*c*a2)+(w'2*Ml*M2*e*c*f))/((a2*h+a3*l)/vv^2- 
(Ml*d*r2+M2*r2*h)); 

g5=((a2*h+a3*l)/w^2+(M2*e*f*h))/((a2*h+a3*l)/w'2-(Ml;i:d*r2+M2*r2*h)); 

g6=(-a2*(Ml*d+M2*h)-(M2*e*a3)+(w>v2*Ml*M2*e"2*d))/((a2*h+a3*l)/w^2- 

(Ml*d*r2+M2*f2*h)); 
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t(l,l)=gl+g3*((g4+gl*g6)/(l-g3*g6)); 
t(l,2)=g2+g3*((g5+g2*g6)/(l-g3*g6)); 
t(2,l)=(g4+gl*g6)/(l-g3*g6); 
t(2,2)=(g5+g2*g6)/(l-g3*g6); 
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%       STRING2AL 

%   This  file  calculates  the  transfer  matrix  for  a  string  of 

length  L, 
%    located  on  the  left  side  of  the  absorber 

%Initial  data 

T=5000; 

L=(7/16) *50; 

P=1.0; 

c=sqrt (T/p) ; 

k=w/c; 

t(l,l)=cos(k*L) ; 

t(l,2)=sin(k*L) / (k*T) ; 

t (2, l)=-k*T*sin(k*L) ; 

t (2,2)=cos (k*L) ; 
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%  STRING2AR 


%  This  file  calculates  the  transfer  matrix  for  a  string  of  length  L, 
%    located  on  the  right  side  of  the  absorber. 

%Initial  data 

T=5000; 

L=(9/16)*50; 

P=1.0; 

c=sqrt(T/p); 

k=w/c; 

t(l,l)=cos(k*L); 

t(l,2)=sin(k*L)/(k*T); 

t(2,l)=-k*T*sin(k*L); 

t(2,2)=cos(k*L); 
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%  STRING4A 

%  This  file  plots  the  overall  T(l,2)  for  a  system  with  3  components. 

W=0; 

Ttotl2=0; 

R=0; 

for  n=  1:200; 

w=.l*n; 

W(n)=w; 

string2al 

Tl=t; 

string2ar 

T3=t; 

string  la 

T2=t; 

Ttot=Tl*T2*T3; 

Ttotl2(n)=Ttot(l,2); 

end 

plot(W,abs(Ttotl2)) 

grid 

xlabel('w(rad/sec)') 

ylabel('T(l,2)') 
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%  STR1NG6A 

%         THIS  FILE  PLOTS  THE  MODE  SHAPES  OF  A  STRING 

%         WITH  AN  OFF  CENTER  DAMPER.  ENTER  FREQ. 

%         FROM  THE  MATLAB  WINDOW  AND  DAMPING 

%         IN  STRING  1 A 

x=0; 

y=0; 

T=5000; 

p=1.0; 

c=sqrt(T/p); 

k=w/c; 

x=0:.  1:50.8; 

for  n=  1:220 

tl(l,l)=cos(k*x(n)); 

tl(l,2)=sin(k*x(n))/(k*T); 

tl(2,l)=-k*T*sin(k*x(n)); 

tl(2,2)=cos(k*x(n)); 

y(n)=tl(l,2); 

f(n)=tl(2,2); 

end 


string  la 

ts(l,l)=cos(21.9*k); 

ts(l,2)=sin(21.9*k)/(k*T); 

ts(2,2)=cos(21.9*k); 

ts(2,l)=-k*T*sin(21.9*k); 

for  n=227:509 

tr(l,l)=cos(k*(x(n)-22.7)); 

tr(  1 ,2)=sin(k*(x(n)-22.7))/(k*T); 

tr(2,l)=-k*T*sin(k*(x(n)-22.7)); 

tr(2,2)=cos(k*(x(n)-22.7)); 

tf=ts*t*tr; 

y(n)=tf(l,2); 

f(n)=tf(2,2); 

end 

plot(x(  l:220),real(y(  1 :220)),x(227:509),real(y(227:509))) 

grid 

xlabelfx  location') 
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ylabel('real(displacement  mode  shape)') 

title('4TH  MODE  MODE  SHAPE,  DAMPER  AT  3L/8') 
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%  STRING  14A 

%  This  file  computes  T(l,2)  for  a  cable  with  an  off  center  absorber. 

function  T12=stringl4a(w) 

string2al 

Tl=t; 

string2ar 

T3=t; 

string  la 

T2=t; 

Ttot=Tl*T2*T3; 

T12=Ttot(l,2); 
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%  MULLER 

function  [zros]=muller(fn,npre v,np,maxit,ep  1  ,ep2,fnreal,zros); 
%  zeros=muller(fn,nprev,maxit,epl,ep2,fnreal,zeros) 
%  This  function  is  called  by  the  driver  program  mulrun.m 
%  Determines  up  to  np  zeros  of  the  function  specified  by  fn, 
%  using  quadratic  interpolation,  i.e.,  Muller's  Method. 
%  This  function  is  a  translation  of  the  Fortran  routine  given 
%  by  S.  D.  Conte  and  Carl  de  Boor,  "Elementary  Numerical 
%  Analysis",  McGraw  Hill  Book  Company,  1980.  The  conversion 
%  was  made  by  Howard  Wilson  and  Chris  Roberts,  Engineering 
%  Mechanics  Department,  University  of  Alabama. 
% 

%  To  find  more  than  one  zero,  function  muller  uses  a  procedure 
%  known  as  deflation,  which  is  implemented  in  function  dflate. 
%  fn:  fuction  that  we  want  to  find  the  roots. 

%  zros:  an  array  of  initial  estimates  of  all  desired  roots,  set  to  zero  if 
%       no  better  estimates  are  available. 
%  epl:  relative  error  tolerance  on  the  root. 
%  ep2:  error  tolerance  on  the  function. 
%  maxit:  maximum  number  of  iterations  per  root  allowed. 
%  np:  number  of  roots  desired. 

%  nprev:  number  of  roots  previously  computed,  normally  zero. 
%  fnreal:  a  logical  number;  if  TRUE,  the  program  forces  all  approximations  to 
%  all  the  roots  to  be  real.  This  makes  it  possible  to  use  this  routine 

%  even  if  f(x)  is  defined  only  for  real  x. 

%  Here  we  can  take  fnreal=l,  if  we  only  want  to  fine  the  real  roots;  we 

%  can  take  fnreal=0  if  we  want  to  find  both  the  real  and  complex  roots. 

fn='stringl4a'; 

eps  1 =max(ep  1 , 1  .e- 1 2);  eps2=max(ep2, 1  .e-20); 
for  i=nprev+l:np;  kount=0; 
while  1 

%    Compute  the  first  three  estimates  for  zero  as 
%    zero,  zero+0.5,  zero-0.5 

zero=zros(i);  h=0.00001+j*0.00001;  zp=zero+h;  zm=zero-h; 

[fzr,dvdflp,kount,zros]=dflate(fn,zp,i,kount,zros); 

[fzr,fzrprv,kount,zros]=dflate(fn,zm,i,kount,zros); 

hprev=-2.0*h;  dvdflp=(fzrprv-dvdflp)/hprev; 

[fzr,fzrdfl,kount,zros]=dflate(fn,zero,i,kount,zros); 

while  2 

%      Do  while  kount<maxit  or  h  is  relatively  big  or  fzr=f(zero) 
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%      is  not  small  or  fzrdfl=fdeflated(zero)  is  not  small  or  not 
%      much  bigger  than  its  previous  value  fzrprv. 

di vdf  1  =(fzrdfl-fzrprv)/h;  di vdf2=(divdf  1  -dvdf  1  p)/(h+hpre v) ; 
hprev=h;  dvdf  lp=di vdf  1;  c=divdfl+h*divdf2; 
sqr=c*c-4.0*fzrdfl*divdf2; 

if  fnreal*real(sqr)<0.0,  sqr=0.0;  else,  sqr=sqrt(sqr);  end; 
if  real(c)*real(sqr)+imag(c)*imag(sqr)<0,  den=c-sqr; 
else,  den=c+sqr;  end 
if  abs(den)<0.0,  den=1.0;  end 
h=-2.0*fzrdfi7den;  fzrprv=fzrdfl;  zero=zero+h; 
if  kount>maxit,  break,  else 
while  3 
[fzr,fzrdfl,kount,zros]=dfiate(fn,zero,i,kount,zros); 
if  zros(i)==zero*  1.001,  break;  end 
%  Check  for  convergence 

if  abs(h)<epsl*abs(zero),  break;  end 
if  max(abs(fzr),abs(fzrdfl))<eps2,  break;  end 
%  Check  for  divergence 

if  abs(fzrdfl)>10.0*abs(fzrprv),  h=h/2;  zero=zero-h; 
else,  break;  end 
end;  %  while  3 
end;  %  end  kount  <=  maxit 
if  zros(i)==zero*  1.001,  break;  end 
if  abs(h)<epsl*abs(zero),  break;  end 
if  max(abs(fzr),abs(fzrdfi))<eps2,  break;  end 
end;  %  while  2 
if  kount>maxit,  break;  end 
if  abs(h)<epsl:f:abs(zero),  break;  end 
if  max(abs(fzr),abs(fzrdfi))<eps2,  break;  end 
end;  %  while  1 
zros(i)=zero; 
end;  %  for  loop 


75 


%  MULRUN 


x=0; 

y=0; 

global  R 

%  input  specs  for  Muller 

forn=l:21 

R=10000*(n-1); 

x(n)=R; 

%    tolerance,  name  of  function  that  computes  T12  etc: 

zrs  =  10.88;  %  Starting  guess  for  the  root,  can  be  real  or  complex 

fn  =  'string  14a';  %  name  of  program  that  computes  T12 
maxit  =  500;  %  max  #  of  iterations 

epl       =       le-5;  %tolerance       on       the       root       (start       with       lower 

%  values,  le-2or  so.) 

tolfcn  =  le-5;  %  relative  tolerance  on  T12 

fnreal  =  0;  %  =  0  complex  search,  =  1,  real  search 

nprev  =  0;  %  leave  these  as  0 

np  =  1 ;  %  and  1 


ep2  =  tolfcn*stringl4a(zrs);  %  This  sets  the  tolerance  on  T12  to  be 

%       le-5      of      T12      estimated      at      your      starting 
%  guess 

%  Call  Muller 

root  =  muller(fn,nprev,np,maxit,epl,ep2,fnreal,zrs); 
y(n)=root; 
end 
zeta=imag(y)./real(y); 

subplot(311) 

plot(x,imag(y)) 

grid 

title('4th  mode') 

xlabel('R') 

ylabel('imag(nat  freq)') 

subplot(312) 

plot(x,real(y)) 
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grid 

xlabel('R') 

ylabel('real(nat  freq)') 

subplot(313) 

plot(x,zeta) 

grid 

xlabel('R') 

ylabel(' modal  damping  ratio') 
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%  DFLATE 

%  This  is  the  subroutine  needed  by  muller. 

function  [fzero,fzrdfl,kount,zros]=dflate(F,zero,i,kount,zros); 

%  (*fzero,fzrdfl,kount,zeros*)=dflate(F,zero,i,kount,zeros) 

%  This  function  which  is  called  by  function  muller  is  translated 

%  from  the  Fortran  routine  given  by  S.  D.  Conte  and  Carl  de  Boor, 

%  "Elementary  Numerical  Analysis",  McGraw  Hill  Book  Company,  1980. 

%  The  translation  was  made  by  Howard  Wilson  and  Chris  Roberts, 

%  Engineering  Mechanics  Department,  University  of  Alabama. 

kount=kount+l;  fzero=feval(F,zero);  fzrdfl=fzero; 
if  i<2,  return;  end 
for  j=2:i 

den=zero-zros(j-l); 

if  abs(den)==0.0 
zros(i)=zero*  1.001;  return 

else 
fzrdfl=fzrdfl/den; 

end 
end 
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MULTABS 

The  following  is  a  typical  sequence  for  using  the  model  MULTABS  to  determine 
system  natural  frequencies  and  mode  shapes  and  optimize  system  damping: 

1)  In  MATLAB,  change  directory  to  MULTABS. 

2)  Edit  STRING  1 A  for  desired  absorber  properties  and  STRTNG2AI  and  STRING2AO 
for  desired  cable  properties. 

3)  In  STRING4A  and  STRING1  A,  ensure  damping  constant  is  set  to  zero.  Run 
STRING4A  to  plot  t!2  of  [Toveraii]  versus  frequency.  The  zeroes  of  this  plot  are  the 
undamped  system  natural  frequencies. 

4)  If  desired,  plot  undamped  system  mode  shapes  as  follows: 

a)  For  mode  of  interest,  enter  natural  frequency  in  MATLAB  window. 

b)  Run  STRTNG6A  to  plot  mode  shape. 

c)  Repeat  for  other  modes  of  interest. 

5)  Perform  absorber  optimization  for  individual  modes  as  follows: 

a)  Edit  MULRUN  to  iterate  over  desired  range  of  damping  constants. 

b)  Also  in  MULRUN,  set  zrs  to  the  undamped  natural  frequency  for  the  mode  of 
interest.  This  serves  as  a  starting  guess  for  the  complex  root  solver. 

6)  Set  damping  constant  to  global  in  STRING  1  A.  Run  MULRUN  to  plot  real  and 
imaginary  natural  frequency  and  damping  ratio  versus  damping  constant.  The  optimum 
damping  constant  is  that  which  maximizes  damping  ratio. 

7)  Plot  damped  system  mode  shape  as  follows: 

a)  Set  damping  constant  to  the  optimum  value  in  STRING1  A. 

b)  Enter  the  complex  natural  frequency  in  the  MATLAB  window. 

c)  Run  STRTNG6A  to  plot  mode  shape. 

8)  Repeat  steps  5  through  7  for  other  modes  of  interest. 


79 


%  STRING  1 A 

%  This  file  contains  the  transfer  matrix  for  a  vibration 
%  absorber. 

%  Initial  data 

T=5000; 

a=.03; 

b=.03; 

c=.2; 

d=.2; 

e=.2; 

f=.2; 

Ml=9; 

M2=M1; 

H=5000; 

11=  12; 

12=11; 

K=10; 

global  R; 

h=c+d; 

l=e+f; 

al=(T*d)+(H*c)+(K*a~2)-(w~2*Il)+(i*w*R*b~2); 

a2=(T*0+(H*e)+(K*a~2)-(w~2*I2Mi*w*R*b~2); 
a3=(K*a~2)+(i*w*R*b~2); 

gl=(-(Ml*l+M2*f)*al- 

(Ml*c*a3)+(w^2*Ml*M2*f*c'2))/((M2*e*al)+(Ml*d+M2*h)*a3+(w'2*Ml*M2*c*e*d)); 

g2=(-(M2*f*hA2)- 
(Ml*l*d~2)+(al*l+a3*h)/w^2)/((M2*e*al)+(Ml*d+M2*h)*a3+(wA2*Ml*M2*c*e*d)); 

g3=- 
((Ml*l*d*c)+(al*l+a3*h)/w^2)/((M2*e*al)+(Ml*d+M2*h)*a3+(w^2*Ml*M2*c*e*d)); 

g4==-((Ml*l+M2*f)*a3+(Ml*c*a2)+(w"2*Ml*M2*e*c*f))/((a2*h+a3*l)/wA2- 
(Ml*d*r2+M2*r2*h)); 

g5=((a2*h+a3*l)/w*2+(M2*e*f*h))/((a2*h+a3*l)/w*2-(Ml*d*r2+M2*r2*h)); 

g6=(-a2*(Ml*d+M2*h)-(M2*e*a3)+(w^2*Ml*M2*e^2*d))/((a2*h+a3*l)/w^2- 
(Ml*d*r2+M2*r2*h)); 


so 


t(l,l)=gl+g3*((g4+gl*g6)/(l-g3*g6)); 
t(l,2)=g2+g3*((g5+g2*g6)/(l-g3*g6)); 
t(2,l)=(g4+gl*g6)/(l-g3*g6); 
t(2,2)=(g5+g2*g6)/(l-g3*g6); 


si 


%  STRING2AI 

%  This  file  calculates  the  transfer  matrix  for  a  string  of  length  L/3, 
%    located  on  the  left  side  of  the  abosorber 

%Initial  data 
T=5000; 
L=(l/3)*50; 
p=1.0; 

c=sqrt(T/p); 

k=w/c; 

t(l,l)=cos(k*L); 

t(l,2)=sin(k*L)/(k*T); 
t(2,l)=-k*T*sin(k*L); 
t(2,2)=cos(k*L); 
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%  STRING2A0 


%  This  file  calculates  the  transfer  matrix  for  a  string  of  length  L/6, 
%    located  on  the  left  side  of  the  abosorber 


%Initial  data 

T=5000; 

L=(l/6)*50; 

p=1.0; 

c=sqrt(T/p); 

k=w/c; 

t(l,l)=cos(k*L); 

t(l,2)=sin(k*L)/(k*T); 

t(2,l)=-k*T*sin(k*L); 

t(2,2)=cos(k*L); 
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%  STRING4A 

%  This  file  plots  the  overall  T(l,2)  for  a  system  with  absorbers  located 
%  at  the  mode  3  antinodes.  Before  running,  set  R  to  the  desired  value 
%  in  STRING  1  A. 

W=0; 
Ttotl2=0; 

R=0; 

for  n=  1:100; 

w=100+.l*n; 

W(n)=w; 

string2ao 

Tl=t; 

string2ai 

T3=t; 

string  la 

T2=t; 

Ttot=Tl  *T^*T3*T^*T3*T2*T1  ° 

Ttotl2(n)=Ttot(l,2); 

end 

plot(W,abs(Ttotl2)) 

grid 

xlabel( '  w(rad/sec) ' ) 

ylabel('T(l,2)') 
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%  STRING6A 

%         THIS  FILE  PLOTS  THE  MODE  SHAPES  OF  A  STRING 

%         WITH  ABSORBERS  AT  THE  MODE  3  ANTINODES. 

%  ENTER  FREQUENCY 

%         FROM  THE  MATLAB  WINDOW  AND  DAMPING 

%  IN  STRING  1 A 

x=0; 

y=0; 

T=5000; 

p=L0; 

c=sqrt(T/p); 

k=w/c; 

x=0:.  1:52.4; 

for  n=  1:84 
tl(l,l)=cos(k*x(n)); 
tl(l,2)=sin(k*x(n))/(k*T); 
tl(2,l)=-k*T*sin(k*x(n)); 

tl(2,2)=cos(k*x(n)); 
y(n)=tl(l,2); 
f(n)=tl(2,2); 
end 


string  la 


tso(l,l)=cos(8.33*k); 

tso(  1 ,2)=sin(8.33*k)/(k*T); 

tso(2,2)=cos(8.33*k); 

tso(2,l)=-k*T*sin(8.33*k); 

tsi(l,l)=cos(16.67*k); 

tsi(l,2)=sin(16.67*k)/(k*T); 

tsi(2,2)=cos(16.67*k); 

tsi(2,l)=-k*T*sin(16.67*k); 

for  n=92:259 

tlc(l,l)=cos(k*(x(n)-9.1)); 

tlc(  1 ,2)=sin(k*(x(n)-9. 1  ))/(k*T); 

tlc(2,l)=-k*T*sin(k*(x(n)-9.1)); 

tlc(2,2)=cos(k*(x(n)-9.1)); 

tf=tso*t*tlc; 

y(n)=tf(l,2); 
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f(n)=tf(2,2); 
end 


for  n=267:434 

trc(l,l)=cos(k*(x(n)-26.6)); 

trc(l,2)=sin(k*(x(n)-26.6))/(k*T); 

trc(2,l)=-k*T*sin(k*(x(n)-26.6)); 

trc(2,2)=cos(k*(x(n)-26.6)); 

tf=tso*t*tsi*t*trc; 

y(n)=tf(l,2); 

f(n)=tf(2,2); 

end 

for  n=442:525 

tr(l,l)=cos(k*(x(n)-44.1)); 

tr(l,2)=sin(k*(x(n)-44.1))/(k*T); 

tr(2,l)=-k*T*sin(k*(x(n)-44.1)); 

tr(2,2)=cos(k*(x(n)-44. 1)); 

tf=tso*t*tsi*t*tsi*t*tr; 

y(n)=tf(l,2); 

f(n)=tf(2,2); 

end 

plot(x(l:84),real(y(l:84)),x(92:259),real(y(92:259)),x(267:434),real(y(267:434)),x(442:525),real(y(442:5 

grid 

xlabel('x  location') 

ylabel('real(displacement  mode  shape)') 

title('3rd  MODE  MODE  SHAPE,  DAMPERS  AT  ANTINODES') 
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%  STRING  14A 

%  This  file  computes  T(l,2)  for  a  cable  with  absorbers  located  at  the 
%    mode  3  antinodes. 

function  T12=stringl4a(w) 

string2ao 

Tl=t; 

string2ai 

T3=t; 

string  la 

T2=t; 

Xtot=Tl*T^*T3*T2*T3*T^*Tl' 
T12=Ttot(l,2); 
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%  MULLER 

function  [zros]=muller(fn,npre v,np,maxit,ep  1  ,ep2,fnreal,zros); 
%  zeros=muller(fn,nprev,maxit,ep  1  ,ep2,fnreal,zeros) 
%  This  function  is  called  by  the  driver  program  mulrun.m 
%  Determines  up  to  np  zeros  of  the  function  specified  by  fn, 
%  using  quadratic  interpolation,  i.e.,  Muller's  Method. 
%  This  function  is  a  translation  of  the  Fortran  routine  given 
%  by  S.  D.  Conte  and  Carl  de  Boor,  "Elementary  Numerical 
%  Analysis",  McGraw  Hill  Book  Company,  1980.  The  conversion 
%  was  made  by  Howard  Wilson  and  Chris  Roberts,  Engineering 
%  Mechanics  Department,  University  of  Alabama. 
% 

%  To  find  more  than  one  zero,  function  muller  uses  a  procedure 
%  known  as  deflation,  which  is  implemented  in  function  dflate. 
%  fn:  fuction  that  we  want  to  find  the  roots. 

%  zros:  an  array  of  initial  estimates  of  all  desired  roots,  set  to  zero  if 
%        no  better  estimates  are  available. 
%  epl:  relative  error  tolerance  on  the  root. 
%  ep2:  error  tolerance  on  the  function. 
%  maxit:  maximum  number  of  iterations  per  root  allowed. 
%  np:  number  of  roots  desired. 

%  nprev:  number  of  roots  previously  computed,  normally  zero. 
%  fnreal:  a  logical  number;  if  TRUE,  the  program  forces  all  approximations  to 
%  all  the  roots  to  be  real.  This  makes  it  possible  to  use  this  routine 

%  even  if  f(x)  is  defined  only  for  real  x. 

%  Here  we  can  take  fnreal=l,  if  we  only  want  to  fine  the  real  roots;  we 

%  can  take  fnreal=0  if  we  want  to  find  both  the  real  and  complex  roots. 

fn=' string  14a'; 

eps l=max(ep  1, l.e-12);  eps2=max(ep2, 1  .e-20); 
for  i=nprev+l:np;  kount=0; 
while  1 

%    Compute  the  first  three  estimates  for  zero  as 
%    zero,  zero+0.5,  zero-0.5 

zero=zros(i);  h=0.00001+j*0.00001;  zp=zero+h;  zm=zero-h; 

[fzr,dvdflp,kount,zros]=dfiate(fn,zp,i,kount,zros); 

[fzr,fzrprv,kount,zros]=dflate(fn,zm,i,kount,zros); 

hprev=-2.0*h;  dvdflp=(fzrprv-dvdflp)/hprev; 

[fzr,fzrdfl,kount,zros]=dflate(fn,zero,i,kount,zros); 

while  2 

%      Do  while  kount<maxit  or  h  is  relatively  big  or  fzr=f(zero) 
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%      is  not  small  or  fzrdfl=fdeflated(zero)  is  not  small  or  not 
%      much  bigger  than  its  previous  value  fzrprv. 

divdfl=(fzrdfl-fzrprv)/h;divdf2=(divdfl-dvdflp)/(h+hprev); 
hprev=h;  dvdflp=divdfl;  c=divdfl+h*divdf2; 
sqr=c*c-4.0*fzrdfl*divdf2; 

if  fnreal*real(sqr)<0.0,  sqr=0.0;  else,  sqr=sqrt(sqr);  end; 
if  real(c)*real(sqr)+imag(c)*imag(sqr)<0,  den=c-sqr; 
else,  den=c+sqr;  end 
if  abs(den)<0.0,  den=1.0;  end 
h=-2.0*fzrdfl/den;  fzrprv=fzrdfl;  zero=zero+h; 
if  kount>maxit,  break,  else 
while  3 
[fzr,fzrdfl,kount,zros]=dflate(fn,zero,i,kount,zros); 
if  zros(i)==zero*  1.001,  break;  end 
%  Check  for  convergence 

if  abs(h)<epsl*abs(zero),  break;  end 
if  max(abs(fzr),abs(fzrdfi))<eps2,  break;  end 
%  Check  for  divergence 

if  abs(fzrdfi)>10.0*abs(fzrprv),  h=h/2;  zero=zero-h; 
else,  break;  end 
end;  %  while  3 
end;  %  end  kount  <=  maxit 
if  zros(i)==zero*  1.001,  break;  end 
if  abs(h)<epsl*abs(zero),  break;  end 
if  max(abs(fzr),abs(fzrdfi))<eps2,  break;  end 
end;  %  while  2 
if  kount>maxit,  break;  end 
if  abs(h)<epsl*abs(zero),  break;  end 
if  max(abs(fzr),abs(fzrdfi))<eps2,  break;  end 
end;  %  while  1 
zros(i)=zero; 
end;  %  for  loop 
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%  MULRUN 


x=0; 

y=0; 

global  R 

%  input  specs  for  Muller 

for  n=  1:21 

R=500*(n-1); 

x(n)=R; 

%    tolerance,  name  of  function  that  computes  T12  etc: 

zrs  =  104;  %  Starting  guess  for  the  root,  can  be  real  or  complex 

fn  =  'string  14a';  %  name  of  program  that  computes  T12 
maxit  =  500;  %  max  #  of  iterations 

epl       =       le-5;  %tolerance       on       the       root       (start       with       lower 

%  values,  le-2or  so.) 

tolfcn  =  le-5;  %  relative  tolerance  on  T12 

fnreal  =  0;  %  =  0  complex  search,  =  1 ,  real  search 

nprev  =  0;  %  leave  these  as  0 

np  =  1 ;  %  and  1 


ep2  =  tolfcn*stringl4a(zrs);  %  This  sets  the  tolerance  on  T12  to  be 

%       le-5      of      T12      estimated      at      your      starting 
%  guess 

%  Call  Muller 

root  =  muller(fn,nprev,np,maxit,epl,ep2,fnreal,zrs); 
y(n)=root; 
end 
zeta=imag(y)./real(y); 

subplot(311) 

plot(x,imag(y)) 

grid 

titlef  mode,  dampers  at  mode  3  antinodes') 

xlabel('R') 

ylabel('imag(nat  freq)') 

subplot(312) 

plot(x,real(y)) 
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grid 

xlabel('R') 

ylabel('real(nat  freq)') 

subplot(313) 

plot(x,zeta) 

grid 

xlabel('R') 

ylabel(' modal  damping  ratio') 


%  DFLATE 

%  This  is  the  subroutine  needed  by  muller. 

function  [fzero,fzrdfl,kount,zros]=dflate(F,zero,i,kount,zros); 

%  (*fzero,fzrdfl,kount,zeros*)=dflate(F,zero,i,kount,zeros) 

%  This  function  which  is  called  by  function  muller  is  translated 

%  from  the  Fortran  routine  given  by  S.  D.  Conte  and  Carl  de  Boor, 

%  "Elementary  Numerical  Analysis",  McGraw  Hill  Book  Company,  1980. 

%  The  translation  was  made  by  Howard  Wilson  and  Chris  Roberts, 

%  Engineering  Mechanics  Department,  University  of  Alabama. 

kount=kount+l;  fzero=feval(F,zero);  fzrdfi=fzero; 
if  i<2,  return;  end 
for  j=2:i 

den=zero-zros(j-l); 

if  abs(den)==0.0 
zros(i)=zero*  1.001;  return 

else 
fzrdfi=fzrdfi/den; 

end 
end 
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HANG 

The  following  is  a  typical  sequence  for  using  the  model  HANG  to  determine 
system  natural  frequencies  and  mode  shapes  and  optimize  system  damping: 

1)  In  MATLAB,  change  directory  to  HANG. 

2)  Edit  STRING1  for  desired  absorber  properties  and  STRING2  for  desired  cable 
properties. 

3)  In  STRING4  and  STRING1,  ensure  damping  constant  is  set  to  zero.  Run  STRING4 
to  plot  t^of  [T0veraii]  versus  frequency.  The  zeroes  of  this  plot  are  the  undamped  system 
natural  frequencies. 

4)  If  desired,  plot  undamped  system  mode  shapes  as  follows: 

a)  For  mode  of  interest,  enter  natural  frequency  in  MATLAB  window. 

b)  Run  STRTNG6  to  plot  mode  shape. 

c)  Repeat  for  other  modes  of  interest. 

5)  Perform  absorber  optimization  for  individual  modes  as  follows: 

a)  Edit  MULRUN  to  iterate  over  desired  range  of  damping  constants. 

b)  Also  in  MULRUN,  set  zrs  to  the  undamped  natural  frequency  for  the  mode  of 
interest.  This  serves  as  a  starting  guess  for  the  complex  root  solver. 

6)  Set  damping  constant  to  global  in  STRING  1.  Run  MULRUN  to  plot  real  and 
imaginary  natural  frequency  and  damping  ratio  versus  damping  constant.  The  optimum 
damping  constant  is  that  which  maximizes  damping  ratio. 

7)  Plot  damped  system  mode  shape  as  follows: 

a)  Set  damping  constant  to  the  optimum  value  in  STRING  1. 

b)  Enter  the  complex  natural  frequency  in  the  MATLAB  window. 

c)  Run  STRING6  to  plot  mode  shape. 

8)  Repeat  steps  5  through  7  for  other  modes  of  interest. 


%  STRING  1 

%  This  file  contains  the  transfer  matrix  for  a  hanging  mass  spring 
%  dashpot  vibration  absorber. 

%  Initial  data 

M=5; 

K=68.45; 

R=15; 


t(l,D=l; 

t(l,2)=0; 

t(2J)=-(w~2*M*(i*w*R+K))/(K-w~2*M+i*w*R); 

t(2,2)=l; 
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%  STRING2 

%  This  file  calculates  the  transfer  matrix  for  a  string  of  length  L. 

%Initial  data 

T=5000; 

L=25; 

p=1.0; 

c=sqrt(T/p); 

k=w/c; 

t(l,l)=cos(k*L); 

t(l,2)=sin(k*L)/(k*T); 
t(2,l)=-k*T*sin(k*L); 
t(2,2)=cos(k*L); 
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%  STRING4 

%  This  file  plots  the  overall  T(l,2)  for  a  system  with  3  components. 

W=0; 

Ttotl2=0; 

R=0; 

for  n=  1:200; 

w=0.1*n; 

W(n)=w; 

string2 

Tl=t; 

T3=t; 

string  1 

T2=t; 

Ttot=Tl*T2*T3; 

Ttotl2(n)=Ttot(l,2); 

end 

plot(W,abs(Ttotl2)) 

grid 

xlabel('w(rad/sec)') 

ylabel('T(l,2)') 
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%  STRING6 

%         THIS  FILE  PLOTS  THE  MODE  SHAPES  OF  A  STRING 

%         WITH  A  MASS  SPRING  DASHPOT  DAMPER  AT  ITS  %  CENTER. 

BEFORE  RUNNING, 

%         ENTER  FREQUENCY  IN  THE  MATLAB  WINDOW  AND 

%         ENTER  DAMPING  IN  STRING  1 . 

y=0; 

f=0; 

T=5000; 

p=1.0; 

c=sqrt(T/p); 

k=w/c; 

x=0:.l:50; 

for  n=  1:251 

tl(l,l)=cos(k*x(n)); 

tl(l,2)=sin(k*x(n))/(k*T); 

tl(2,l)=-k*T*sin(k*x(n)); 

tl(2,2)=cos(k*x(n)); 

y(n)=tl(l,2); 

f(n)=tl(2,2); 

end 


string  1 

ts(l,l)=cos(25*k); 

ts(l,2)=sin(25*k)/(k*T); 

ts(2,2)=cos(25*k); 
ts(2,l)=-k*T*sin(25*k); 

forn=252:501 

tr(l,l)=cos(k*(x(n)-25)); 

tr(  1 ,2)=sin(k*(x(n)-25))/(k*T); 

tr(2,l)=-k*T*sin(k*(x(n)-25)); 

tr(2,2)=cos(k*(x(n)-25)); 

tf=ts*t*tr; 

y(n)=tf(l,2); 

f(n)=tf(2,2); 

end 

plot(x(l:251),real(y(l:251)),x(252:501),real(y(252:501))) 

grid 
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xlabelfx  location') 

ylabel('real(displacement  mode  shape)') 
title('3rd  MODE  MODE  SHAPE') 
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%  STRING  14 

% 


%  This  file  computes  T(l,2)  for  a  system  with  3  components. 


function  T12=stringl4(w) 

string2 

Tl=t; 

T3=t; 

string  1 

T2=t; 

Ttot=Tl*T2*T3; 

T12=Ttot(l,2); 
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%  MULLER 

function  [zros]=muller(fn,nprev,np,maxit,epl,ep2,fnreal,zros); 
%  zeros=muller(fn,nprev,maxit,epl,ep2,fnreal,zeros) 
%  This  function  is  called  by  the  driver  program  mulrun.m 
%  Determines  up  to  np  zeros  of  the  function  specified  by  fn, 
%  using  quadratic  interpolation,  i.e.,  Muller's  Method. 
%  This  function  is  a  translation  of  the  Fortran  routine  given 
%  by  S.  D.  Conte  and  Carl  de  Boor,  "Elementary  Numerical 
%  Analysis",  McGraw  Hill  Book  Company,  1980.  The  conversion 
%  was  made  by  Howard  Wilson  and  Chris  Roberts,  Engineering 
%  Mechanics  Department,  University  of  Alabama. 
% 

%  To  find  more  than  one  zero,  function  muller  uses  a  procedure 
%  known  as  deflation,  which  is  implemented  in  function  dflate. 
%  fn:  fuction  that  we  want  to  find  the  roots. 

%  zros:  an  array  of  initial  estimates  of  all  desired  roots,  set  to  zero  if 
%        no  better  estimates  are  available. 
%  epl:  relative  error  tolerance  on  the  root. 
%  ep2:  error  tolerance  on  the  function. 
%  maxit:  maximum  number  of  iterations  per  root  allowed. 
%  np:  number  of  roots  desired. 

%  nprev:  number  of  roots  previously  computed,  normally  zero. 
%  fnreal:  a  logical  number;  if  TRUE,  the  program  forces  all  approximations  to 
%  all  the  roots  to  be  real.  This  makes  it  possible  to  use  this  routine 

%  even  if  f(x)  is  defined  only  for  real  x. 

%  Here  we  can  take  fnreal=l,  if  we  only  want  to  fine  the  real  roots;  we 

%  can  take  fnreal=0  if  we  want  to  find  both  the  real  and  complex  roots. 

fn='stringl4'; 

eps l=max(ep  1 , 1  .e- 12);  eps2=max(ep2, 1  .e-20); 
for  i=nprev+l:np;  kount=0; 
while  1 

%    Compute  the  first  three  estimates  for  zero  as 
%    zero,  zero+0.5,  zero-0.5 

zero=zros(i);  h=0.00001+j*0.00001;  zp=zero+h;  zm=zero-h; 

[fzr,dvdflp,kount,zros]=dflate(fn,zp,i,kount,zros); 

[fzr,fzrprv,kount,zros]=dflate(fn,zm,i,kount,zros); 

hprev=-2.0*h;  dvdflp=(fzrprv-dvdflp)/hprev; 

[fzr,fzrdfl,kount,zrosl=dflate(fn,zero,i,kount,zros); 

while  2 

%      Do  while  kount<maxit  or  h  is  relatively  big  or  fzr=f(zero) 


loo 


-2- 


%      is  not  small  or  fzrdfl=fdeflated(zero)  is  not  small  or  not 
%      much  bigger  than  its  previous  value  fzrprv. 

di vdf  1 =(fzrdfl-fzrprv)/h ;  di vdf 2=(di vdf  1  -dvdf  1  p)/(h+hpre v) ; 
hprev=h;  dvdflp=divdfl;  c=divdfl+h*divdf2; 
sqr=c*c-4.0*fzrdfl*divdf2; 

if  fnreal*real(sqr)<0.0,  sqr=0.0;  else,  sqr=sqrt(sqr);  end; 
if  real(c)*real(sqr)+imag(c)*imag(sqr)<0,  den=c-sqr; 
else,  den=c+sqr;  end 
if  abs(den)<0.0,  den=1.0;  end 
h=-2.0*fzrdfi/den;  fzrprv=fzrdfl;  zero=zero+h; 
if  kount>maxit,  break,  else 
while  3 
[fzr,fzrdfl,kount,zros]=dflate(fn,zero,i,kount,zros); 
if  zros(i)==zero*  1.001,  break;  end 
%  Check  for  convergence 

if  abs(h)<epsl*abs(zero),  break;  end 
if  max(abs(fzr),abs(fzrdfl))<eps2,  break;  end 
%  Check  for  divergence 

if  abs(fzrdfi)>10.0*abs(fzrprv),  h=h/2;  zero=zero-h; 
else,  break;  end 
end;  %  while  3 
end;  %  end  kount  <=  maxit 
if  zros(i)==zero*  1.001,  break;  end 
if  abs(h)<epsl*abs(zero),  break;  end 
if  max(abs(fzr),abs(fzrdfi))<eps2,  break;  end 
end;  %  while  2 
if  kount>maxit,  break;  end 
if  abs(h)<epsl*abs(zero),  break;  end 
if  max(abs(fzr),abs(fzrdfl))<eps2,  break;  end 
end;  %  while  1 
zros(i)=zero; 
end;  %  for  loop 
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%  MULRUN 


CI     IK******:!:******:):****:):****************** 

x=0; 

y=0; 

global  R 

%  input  specs  for  Muller 

for  n=  1:41 

R=.5*n; 

x(n)=R; 

%    tolerance,  name  of  function  that  computes  T12  etc: 

zrs  =  4.07;  %  Starting  guess  for  the  root,  can  be  real  or  complex 

fn  =  'string  14';  %  name  of  program  that  computes  T12 
maxit  =  500;  %  max  #  of  iterations 

epl       =       le-5;  %tolerance       on       the       root       (start       with       lower 

%         values,  le-2or  so.) 

tolfcn  =  le-5;  %  relative  tolerance  on  T12 

fnreal  =  0;  %  =  0  complex  search,  =  1 ,  real  search 

nprev  =  0;  %  leave  these  as  0 

np  =  1 ;  %  and  1 


ep2  =  tolfcn*stringl4(zrs);  %  This  sets  the  tolerance  on  T12  to  be 

%       le-5      of      T12      estimated      at      your      starting 
%  guess 

%  Call  Muller 

root  =  muller(fn,nprev,np,maxit,epl,ep2,fnreal,zrs); 
y(n)=root; 
end 
zeta=imag(y)./real(y); 

subplot(311) 

plot(x,imag(y)) 

grid 

title("4th  mode') 

xlabel('R') 

ylabel('imag(nat  freq)') 

subplot(312) 

plot(x,real(y)) 
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grid 

xlabel('R') 

ylabel('real(nat  freq)') 

subplot(313) 

plot(x,zeta) 

grid 

xlabel('R') 

ylabel(' modal  damping  ratio') 
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%  DFLATE 

%  This  is  the  subroutine  needed  by  muller. 

function  [fzero,fzrdfi,kount,zros]=dfiate(F,zero,i,kount,zros); 

%  (*fzero,fzrdfl,kount,zeros*)=dflate(F,zero,i,kount,zeros) 

%  This  function  which  is  called  by  function  muller  is  translated 

%  from  the  Fortran  routine  given  by  S.  D.  Conte  and  Carl  de  Boor, 

%  "Elementary  Numerical  Analysis",  McGraw  Hill  Book  Company,  1980. 

%  The  translation  was  made  by  Howard  Wilson  and  Chris  Roberts, 

%  Engineering  Mechanics  Department,  University  of  Alabama. 

kount=kount+l;  fzero=feval(F,zero);  fzrdfl=fzero; 
if  i<2,  return;  end 
for  j=2:i 

den=zero-zros(j-l); 

if  abs(den)==0.0 
zros(i)=zero*  1.001;  return 

else 
fzrdfl=fzrdfi/den; 

end 
end 
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Appendix  B:  Transfer  Matrix  Library 

The  following  appendix  is  a  direct  reproduction  of  the  transfer  matrix  library 
derived  by  Li  Li  [5].  It  is  included  here  so  that  this  thesis  can  serve  as  a  stand  alone 
reference  on  transfer  matrix  analysis 
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Appendix  A 


Transfer  Matrix  Library 


Definitions  of  the  State  Vectors  and  Transfer  Matrices:  (see  Section  2.1) 

A.l      Strings 

1)  A  string  with  varying  tension 

The  WKB  solution  to  a  string  with  varying-tension  (see  Fig.  A-l)  is  given  in  [20] 


as 


The  transfer  matrix  can  be  found  from  this  solution  as 


(A.l) 


Fi  = 


*11      ^12 
*21       ^22 


(A.2) 


where: 


tn  =  (I^-^cosW  +  ± — ^— rsinW 


t12  = 


K/%TMTi)i 


sinW 
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i-1 


t 


Figure  A-l:  A  string  with  varying  tension 


t22  =  (Tl-)xcosW  -  ± %—T.sinW 


where  W  =  u  J0C  — ^ — .    If  the  tension  varies  linearly,  then  W 
v^]  and  T(  =  TU  =  T. 


^[y/Z 


2)  A  string  with  constant  tension 


For  a  string  with  constant  tension  (see  Fig.  A-2),  with  k  =    yj\p'"    :  — i,  the 
transfer  matrix  is 


Fi  = 


cosklc         -gpsinkl,. 
—  kT 'sinkl,.       cosklc 

where  /c  is  the  distance  between  the  two  state  vectors  at  i  —  1  and  i. 


(A.3) 
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i-1 


T      <- 


/ 


k 


H 


■>     T 


Figure  A-2:  A  string  with  constant  tension 


Ti-i  < 


>    T, 


Figure  A-3:  A  general,  rigid  lump 

A. 2      Rigid  Lumps 

1)  A  general,  rigid  lump 

Section  2.5  gives  the  detailed  derivation  of  the  transfer  matrix  for  a  general,  rigid 
lump  (see  Fig.  A-3),  the  result  is 


Ft  = 


Ti{Le  -  I)  +  ZUJ  -  I^[r2  +  (Lc  -  /)2]  L] 

-Mu2[Ti{Le  -l)  +  Ti.,1  -  MwV]        Ti{Le  - 1)  +  TMl  -  A/u>2[r2  +  I2} 

(A.4) 

where  A  =  Tx{Le  -  I)  +  ^.^  -  A/u;2[r2  -  (Ze  -  /)/]. 
2)  Application  1:  a  uniform  bar 
For  a  rigid,  uniform  bar  (see  Fig.  A-4)  with  length  If,  and  linear  density  pi,  then 


108 


i-1 


\,  * 


/ 


\- 


■>     T 


Figure  A-4:  A  rigid,  uniform  bar 


■>    Ti 


Figure  A-5:  A  uniform,  solid  sphere 


Lc  =  lb,  I  =  -^,  M  =  pihi  and  r2  =  ~.  The  transfer  matrix 


is 


Fi  = 


(A.5) 


4(3  -  /J  l-f 

-Mu>2(12-Ih)    4(3  -A) 

where  /{,  =  iw  and  T  is  the  average  tension  if  there  is  a  difference  in  tension  at  the 
two  ends  of  the  lump. 

3)  Application  2:  a  solid  sphere 

For  a  uniform,  solid  sphere  (see  Fig.  A-5)  with  diameter  D  and  volume  density 
pv,  then  Le  =  D,  I  =  f ,  M  =  =£^=.,  and  r2  =  ££.  The  transfer  matrix  is 


Fi== 


20  +  37, 


20  -  77,  sofi 

-Mw2(20-27,)    20-77, 


(A.6) 
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^ 


*  T; 


^ 


Figure  A-6:  A  cone 


where  /,  =  ^  and  T  is  the  average  tension  if  there  is  a  difference  in  tension  at 
two  ends  of  the  lump. 

4)  Application  3:  a  cone 

For  a  cone  (see  Fig.  A-6)  with  height  h,  bottom  radius  rj,,  and  volume  density  pv, 
then  Lc  =  D,l=  \h,  M  =  ?&£■,  and  r2  =  %&&±.  The  transfer  matrix  is 


Fi  = 


l-/c(r2-/i') 


Ah 


-M^[l-Ic(rl  +  \h')     1-/^+4^) 


(A.7) 


where  Ie  = 


3Mu5 


5)  Application  4:  a  point  mass 

The  transfer  matrix  for  a  point  mass  with  mass  M  can  be  obtained  directly  from 
the  transfer  matrix  for  the  general,  rigid  lump  with  Le  =  /  =  r  =  0: 


&  = 


1        0 

-Mu2    1 


(A.8) 


110 


c  =  l1+l2 


d=l^l4 


e  =  lc*-l 


?l6 


Figure  A-7:  A  general,  rigid  lump  with  damping  devices  (type  I) 

A. 3      Damping  Devices 

For  a  rigid  lump  with  a  damping  device,  the  transfer  matrix  is  defined  the  same 
way  as  that  for  a  rigid  lump  without  a  damping  device, 


y 

Ty' 


^11       ^12 
^21       ^22 


y 

Ty' 


(A.9) 


i-l 


The  expressions  for  the  elements  of  the  transfer  matrix  are  as  follows: 


.               .      74  +  7i7e         .               ,      7s  +  7276        ,  74  +  7i7e 

tn  =7i+73- ,      h2=  72+73- ,       t2i  = 


1  -  737s  1  -  7376 

For  the  system  shown  in  Fig.  A-7,  7's  are  given  as  follows: 


.     _  7s  +  7276 

1  -  7376  '  1  -  7376 


(A.10) 


7i  = 


9s{9i9i2  -  gggn)  -  95J929U  +  9199) 
9i2{9296  -  9a9s)  +  99{g*.gs  +  939e) 


(A.H) 
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ffs(ff9  -  929io)  -  98J9z9io  +  912) 

912(9296  -  g^gs)  +  gsig^gs  +  gzg*) 


(A.12) 


_      99(9z97  -  gs)  +  912(9297  4-  gs)  ,A  13) 

9\2{929&  -  9a9%)  +  99{9a95  +  9s9s) 

9\(9s9\8  -  9sgiz)  +  9n{9295  +  939s)  ,  A  n  ., 

74  = ; ; ; : 7 r        (A. 14) 

9is{9397  ~  9s)  +  g\z{9297  +  9s)  +  9\a\929s  +  53^s) 

08^13  -  9l5{g29s  +  939s)  -  9sgiS /  A    i  r\ 

75  =  } v 7 v -, r  (A. 15) 

g\s(gzg7  -  gs)  +  g\z(g2g7  +  gs)  +  gi4[g2gs  +  gzgs) 
^13(^2^6  -  g*gs)  -  g\s{g2gi  +  ^3^)  +  gisig^gs  +  939e)  ,  A  .  a, 

76  =  -, ; ; ; ; ; —  (A. lb) 

g\s\gzg7  -  gs)  +  ^13(^7  +  g%)  +  g\A{g2gs  +  gzgs) 

where  <7's  are  given  by 

9l  =  u;2M:(l  -  ±),  g2  =  u;2(M^  +  M2lf),  g2  =  ^(M2l-j  +  M3^), 

g4  =  u;2M3(l  -  k),  gh  =  2*.  4.  u;2M3^  +  &   y6  =  ^  +  u,2M3/5(^  -  1), 

g7  =  e,g8  =  %,g9  =  ^+  u;2M^  +  $-, 

=  c,gil  =  ^-u,2MMl-li),9i2  =  §, 

=  J-  u*Mz1^  +  &,  gu  =  U,  gls  =  l3, 

'2M3/4(1  -  lf)  +  %,  gl7  =  &  +  u;2M1/3(l  -  ^),  p18  =  a  -  u;2Mx^  +  & 
*i  =  2i-i/i  +  ffjfe  +  Kia\  -  ^1,  +  juR.bl 
a2  =  Hrl4  +  Hil3  +  K\al  +  A'2a2  +  jw(Rib*  +  £2&2)  -  w2J2 
q3  =  Tt/6  +  Hrls  +  tf2a2  -  w2/3  +  juR2b\ 
ft  =  K\a\  +  juRJl  ft  =  K2a\  +  juR2b\ 

where  Hi  and  Hr  are  the  forces  along  the  structural  length  at  the  left  and  right 
joint,  respectively,  which  can  be  found  by  statics. 

For  the  system  shown  in  Fig.  A-8,  7's  are  given  as  follows: 


-{Mxl  +  M2f)ax  -  Mxcaz  +  u2M1M2fc2  -M2fk2  -  MJd2  +  (aj  +  azh)/u>\ 

71  "     M2eax  +  (Mid  +  M2h)az  +  u2MxM2ced  '      72  "  M7eal  +  (M1d  +  M2h)az  +  uPM-iM** 

(A.17) 


g\o 
g\z 
gie  =  v 
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7s,n  =  ^8^13  +  9s9[3  ~  9is{92g's  +  g^g's)  -  9s9is  -  9s9i8) 

7e,n  =  9[z(9296  -  949s)  +  gizigzg*  -  gM  -  g'isigigs  +  gag*)  -  gxeigig's  +  gzg's) 

+g'is{g<g-o  +  gzge)  +  gisig^g's  +  gM 
ild  =  g'lsig^gr  -g$)-  gisg's  +  ^3(^7  +  g*)  +  gi2g's  +  guigig's  +  gig's) 


ls,d  —  le,d  —  l4,d 


g'5=J^  +  ld)'9*=J"b2-e'9s=J"b7l 

g'9  =  ju,b\\  +  i),  g'u  =  jubn  gi2  =  jujhn 

g[3  =  i^W  +  \),  gis  =  J^2h  g'u  =  M2i, 

g[s=MVi  +  le) 

For  the  system  shown  in  Fig.  A-8,  7's  are  given  by 


,  jub7{[Ml{l+c)+M,f]il,d  +  [Mld+M2(e+h)]il,n}        ,   _   J^{-l^7^[Mld^M3{e+h)]-ri,n} 


'l,d 


,   __   j»V2{i#'/i,d  +  [Mld+M,{e  +  h)]',iin}        ,   __   j^I{-[Mi(:+c)+M2/l7^,i-^74,n} 


73   = 


73.4 


74  = 


7«.j 


75   —  ^3  5  76   _ 


'S.J 


~«,J 


A. 4      Supports 


For  the  support  shown  in  Fig.  A-9,  the  transfer  matrix  is 


Fx  = 


L,-z 


0 


Iquj7  -r,(Lc-z)-r,_iz  z 

(Le-z)z  Lc-2 


,      r^0,Xe 


(A.26) 


where  70  is  the  mass  moment  of  inertia  of  the  rigid  lump  about  the  pinned  point. 
For  the  support  shown  in  Fig.  A-10,  the  transfer  matrix  is 


Ft  = 


Lr-Z 

Z 

loJ2  -T(Lc-z) 


0 


,      z^0,Lc 


(A.27) 


(Le-z)z  Lc-z 

where  I0  is  the  mass  moment  of  inertia  of  the  rigid  lump  about  the  pinned  point. 


Ii: 


Ti-1  < 


>    T: 


Figure  A-9:  An  intermediate  support 


>    T 


Figure  A-10:  An  end  support 
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Appendix  C: 


Analysis  of  a  Single  Degree  of  Freedom  System  Using  the 
Transfer  Matrix  and  Complex  Natural  Frequency  Approach 

The  System: 


►    F, 


*•    x^ 


O  0 

The  Transfer  Matrix  Relates  the  State  Vectors  at  Points  I  and  2: 


Til    Ti: 
T:i    T:: 


Xl 

Fi 


(1) 


Derive  the  Transfer  Matrix  of  the  System  Between  Points  1  and  2: 


-F,  -  K(x2  -  xi)  +  jcoR(x2  -  xi)  =  (K  +jcoR)(x2  -  Xi) 


solving  for  x2: 


x2  =  -F,/(K+joR)  +  x1 

F!+F2=M(d:x2/dr) 

Fi  +  F2  =  -co2Mx2 

writing  F2  in  terms  of  Fi  and  Xi: 


(2) 


F2  =  -F,  -  co-M(-Fi/(K  +  jcoR)  +  Xi) 
F2  =  -co2Mx,  +  (orM/(K  +  jcoR)  -1)F, 


(3) 
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Equations  2  and  3  result  in  the  transfer  matrix: 


Tn   Ti: 
T:i   T:: 


CO 


-l/(K  +  jtf>R) 

:M     -\  +  co:M/(k  +  ]coR) 


Applying  Boundary  Conditions: 

Xi  =  0   due  to  system  being  attached  to  wall 

F2  =  0  because  we  want  to  study  the  free  response  of  the  system 

so,  equation  1  becomes: 


x: 

0 


Tn    Ti: 
T:i    T;: 


0 
Fi 


yields  the  characteristic  equation: 
T22  =  0 

-\  +  co:M/(k  +  jcoR)  =0 

-co2M+jcoR  +  K  =  0 
The  Complex  Frequency  Root: 

applying  quadratic  formula  with  4KM  <  R  : 


co  = 


V4KM-R2      .   R 


2M 


+  J 


2M 


tUnsal   '  Jwiniaginan' 


(4) 


(5) 


From  Section  3.3.1  of  thesis,  a  damped  harmonic  response  can  be  written: 
x(t)  =  Acos(Qreait)e"<o!mag' 


x(t)  =  A  cos 


V4KM-R: 


2M 


exp 


R 

2M 


t       (6) 


16 


Comparison  With  Traditional  Method: 

In  the  traditional  method,  the  characteristic  equation  is  the  same  as 

equation  4: 

-co2M+jc3R  +  K  =  0 

but  the  complex  root  is  not  solved  for.  Instead  the  following  are  defined: 

coo  =  J —  ,  the  undamped  natural  frequency  (7) 

C  =  ,  the  damping  ratio  (8) 

2(OoM 

cod  =  coo  V 1  ~C~    > tne  damped  natural  frequency         (9) 

The  displacement  solution  is  then  written  as: 

x(t)  =  A  cos(  co<it )  exp(-^coot)  (10) 

When  equations  7,  8  and  9  are  used  to  write  equation  10  in  terms  of 

K,  R  and  M.  The  result  is: 


x(t)  =  A  cos 


V4KM-R-    )         (      R 

-t    exp tl     (11) 

2M  HV    2M 


which  is  the  same  as  equation  6,  the  result  obtained  using  the  complex 
frequency  approach. 
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